diff --git a/.conda/meta.yaml b/.conda/meta.yaml
index 7cd02969fe..e781b46049 100644
--- a/.conda/meta.yaml
+++ b/.conda/meta.yaml
@@ -14,7 +14,6 @@ requirements:
- {{ compiler('c') }} # [unix]
host:
- cython >=0.25.2
- - lpsolve55
- numpy
- openbabel >=3
- pydas >=1.0.2
@@ -39,7 +38,6 @@ requirements:
- h5py
- jinja2
- jupyter
- - lpsolve55
- markupsafe
- matplotlib >=1.5
- mopac
diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 4e2cc7a77d..5e389fc77a 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -33,8 +33,8 @@ on:
schedule:
# * is a special character in YAML so you have to quote this string
- cron: "0 8 * * *"
- # runs on all branches on both RMG-Py and forks
- push:
+ # allow running on RMG-Py on a pushed branch, only if triggered manually
+ workflow_dispatch:
# runs on PRs against RMG-Py (and anywhere else, but we add this for RMG-Py)
pull_request:
@@ -51,10 +51,20 @@ env:
jobs:
- build-osx:
- runs-on: macos-13
+ build-and-test:
+ strategy:
+ fail-fast: false
+ matrix:
+ python-version: ["3.9"]
+ os: [macos-13, macos-latest, ubuntu-latest]
+ test_julia: [yes, no]
+ runs-on: ${{ matrix.os }}
+ name: ${{ matrix.os }} Build and Test Python ${{ matrix.python-version }} (Julia? ${{ matrix.test_julia }})
# skip scheduled runs from forks
if: ${{ !( github.repository != 'ReactionMechanismGenerator/RMG-Py' && github.event_name == 'schedule' ) }}
+ env:
+ # This is true only if this is a reference case for the regression testing:
+ REFERENCE_JOB: ${{ github.ref == 'refs/heads/main' && github.repository == 'ReactionMechanismGenerator/RMG-Py' }}
defaults:
run:
shell: bash -l {0}
@@ -62,27 +72,17 @@ jobs:
- name: Checkout RMG-Py
uses: actions/checkout@v4
- # Step to create a custom condarc.yml before setting up conda
- - name: Create custom conda config file
- run: |
- RUNNER_CWD=$(pwd)
- echo "channels:" > $RUNNER_CWD/condarc.yml
- echo " - conda-forge" >> $RUNNER_CWD/condarc.yml
- echo " - rmg" >> $RUNNER_CWD/condarc.yml
- echo " - cantera" >> $RUNNER_CWD/condarc.yml
- echo "show_channel_urls: true" >> $RUNNER_CWD/condarc.yml
-
- # configures the mamba environment manager and builds the environment
- - name: Setup Miniforge Python 3.7
+ - name: Setup Miniforge Python ${{ matrix.python-version }}
uses: conda-incubator/setup-miniconda@v3
with:
environment-file: environment.yml
miniforge-variant: Miniforge3
miniforge-version: latest
- python-version: 3.7
- condarc-file: condarc.yml
+ python-version: ${{ matrix.python-version }}
activate-environment: rmg_env
use-mamba: true
+ show-channel-urls: true
+ channels: conda-forge,cantera,rmg
# list the environment for debugging purposes
- name: mamba info
@@ -109,73 +109,39 @@ jobs:
make clean
make
- build-and-test-linux:
- runs-on: ubuntu-latest
- # skip scheduled runs from forks
- if: ${{ !( github.repository != 'ReactionMechanismGenerator/RMG-Py' && github.event_name == 'schedule' ) }}
- env:
- # This is true only if this is a reference case for the regression testing:
- REFERENCE_JOB: ${{ github.ref == 'refs/heads/main' && github.repository == 'ReactionMechanismGenerator/RMG-Py' }}
- defaults:
- run:
- shell: bash -l {0}
- steps:
- - name: Checkout RMG-Py
- uses: actions/checkout@v4
-
- # Step to create a custom condarc.yml before setting up conda
- - name: Create custom condarc.yml
+ # Setup Juliaup
+ - name: Set Julia paths
+ if: matrix.test_julia == 'yes'
run: |
- RUNNER_CWD=$(pwd)
- echo "channels:" > $RUNNER_CWD/condarc.yml
- echo " - conda-forge" >> $RUNNER_CWD/condarc.yml
- echo " - rmg" >> $RUNNER_CWD/condarc.yml
- echo " - cantera" >> $RUNNER_CWD/condarc.yml
- echo "show_channel_urls: true" >> $RUNNER_CWD/condarc.yml
-
- # configures the mamba environment manager and builds the environment
- - name: Setup Miniforge Python 3.7
- uses: conda-incubator/setup-miniconda@v3
+ # echo "JULIAUP_DEPOT_PATH=$CONDA/envs/rmg_env/.julia" >> $GITHUB_ENV
+ # echo "JULIAUP_DEPOT_PATH=$CONDA/envs/rmg_env/.julia" >> $GITHUB_PATH
+ # echo "JULIA_DEPOT_PATH=$CONDA/envs/rmg_env/.julia" >> $GITHUB_ENV
+ # echo "JULIA_DEPOT_PATH=$CONDA/envs/rmg_env/.julia" >> $GITHUB_PATH
+ # echo "JULIA_CONDAPKG_EXE=$CONDA/condabin/mamba" >> $GITHUB_ENV
+ # echo "JULIA_CONDAPKG_EXE=$CONDA/condabin/mamba" >> $GITHUB_PATH
+ # echo "JULIA_CONDAPKG_EXE=$CONDA/condabin/mamba" >> $GITHUB_ENV
+ # echo "JULIA_CONDAPKG_EXE=$CONDA/condabin/mamba" >> $GITHUB_PATH
+ echo "JULIA_CONDAPKG_BACKEND=Current" >> $GITHUB_ENV
+ # echo "JULIA_CONDAPKG_BACKEND=Current" >> $GITHUB_PATH
+
+ - name: Setup Juliaup
+ if: matrix.test_julia == 'yes'
+ uses: julia-actions/install-juliaup@v2
with:
- environment-file: environment.yml
- miniforge-variant: Miniforge3
- miniforge-version: latest
- python-version: 3.7
- condarc-file: condarc.yml
- activate-environment: rmg_env
- use-mamba: true
+ channel: '1.9'
- # list the environment for debugging purposes
- - name: mamba info
- run: |
- mamba info
- mamba list
-
- # Clone RMG-database
- - name: Clone RMG-database
- run: |
- cd ..
- git clone -b $RMG_DATABASE_BRANCH https://github.com/ReactionMechanismGenerator/RMG-database.git
-
- # modify env variables as directed in the RMG installation instructions
- - name: Set Environment Variables
- run: |
- RUNNER_CWD=$(pwd)
- echo "PYTHONPATH=$RUNNER_CWD/RMG-Py:$PYTHONPATH" >> $GITHUB_ENV
- echo "$RUNNER_CWD/RMG-Py" >> $GITHUB_PATH
-
- # RMG build step
- - name: make RMG
- run: |
- make clean
- make
+ - name: Check Julia version
+ if: matrix.test_julia == 'yes'
+ run: julia --version
# RMS installation and linking to Julia
- name: Install and link Julia dependencies
+ if: matrix.test_julia == 'yes'
timeout-minutes: 120 # this usually takes 20-45 minutes (or hangs for 6+ hours).
+ # JULIA_CONDAPKG_EXE points to the existing conda/mamba to avoid JuliaCall from installing their own. See https://juliapy.github.io/PythonCall.jl/stable/pythoncall/#If-you-already-have-a-Conda-environment.
run: |
- python -c "import julia; julia.install(); import diffeqpy; diffeqpy.install()"
- julia -e 'using Pkg; Pkg.add(PackageSpec(name="ReactionMechanismSimulator",rev="for_rmg")); using ReactionMechanismSimulator'
+ mamba install conda-forge::pyjuliacall
+ julia -e 'using Pkg; Pkg.add(Pkg.PackageSpec(name="ReactionMechanismSimulator", url="https://github.com/hwpang/ReactionMechanismSimulator.jl.git", rev="fix_installation")); using ReactionMechanismSimulator'
- name: Install Q2DTor
run: echo "" | make q2dtor
@@ -192,7 +158,7 @@ jobs:
run: |
for regr_test in aromatics liquid_oxidation nitrogen oxidation sulfur superminimal RMS_constantVIdealGasReactor_superminimal RMS_CSTR_liquid_oxidation RMS_liquidSurface_ch4o2cat fragment RMS_constantVIdealGasReactor_fragment minimal_surface;
do
- if python-jl rmg.py test/regression/"$regr_test"/input.py; then
+ if python rmg.py test/regression/"$regr_test"/input.py; then
echo "$regr_test" "Executed Successfully"
else
echo "$regr_test" "Failed to Execute" | tee -a $GITHUB_STEP_SUMMARY
@@ -282,7 +248,7 @@ jobs:
echo ""
# Compare the edge and core
- if python-jl scripts/checkModels.py \
+ if python scripts/checkModels.py \
"$regr_test-core" \
$REFERENCE/"$regr_test"/chemkin/chem_annotated.inp \
$REFERENCE/"$regr_test"/chemkin/species_dictionary.txt \
@@ -299,7 +265,7 @@ jobs:
cat "$regr_test-core.log" || (echo "Dumping the whole log failed, please download it from GitHub actions. Here are the first 100 lines:" && head -n100 "$regr_test-core.log")
echo " "
echo ""
- if python-jl scripts/checkModels.py \
+ if python scripts/checkModels.py \
"$regr_test-edge" \
$REFERENCE/"$regr_test"/chemkin/chem_edge_annotated.inp \
$REFERENCE/"$regr_test"/chemkin/species_edge_dictionary.txt \
@@ -320,7 +286,7 @@ jobs:
if [ -f test/regression/"$regr_test"/regression_input.py ];
then
echo ""
- if python-jl rmgpy/tools/regression.py \
+ if python rmgpy/tools/regression.py \
test/regression/"$regr_test"/regression_input.py \
$REFERENCE/"$regr_test"/chemkin \
test/regression/"$regr_test"/chemkin
@@ -384,7 +350,7 @@ jobs:
# who knows ¯\_(ツ)_/¯
#
# taken from https://github.com/docker/build-push-action
- needs: build-and-test-linux
+ needs: build-and-test
runs-on: ubuntu-latest
if: github.ref == 'refs/heads/main' && github.repository == 'ReactionMechanismGenerator/RMG-Py'
steps:
diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index ef3509d78f..084224f954 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -26,15 +26,17 @@ jobs:
with:
fetch-depth: 0
- - name: Setup Mambaforge Python 3.7
+ - name: Setup Miniforge Python 3.9
uses: conda-incubator/setup-miniconda@v3
with:
environment-file: environment.yml
miniforge-variant: Miniforge3
miniforge-version: latest
- python-version: 3.7
+ python-version: 3.9
activate-environment: rmg_env
use-mamba: true
+ show-channel-urls: true
+ channels: conda-forge,cantera,rmg
- name: Install sphinx
run: mamba install -y sphinx
@@ -60,12 +62,6 @@ jobs:
make clean
make
- - name: Install and link Julia dependencies
- timeout-minutes: 120 # this usually takes 20-45 minutes (or hangs for 6+ hours).
- run: |
- python -c "import julia; julia.install(); import diffeqpy; diffeqpy.install()"
- julia -e 'using Pkg; Pkg.add(PackageSpec(name="ReactionMechanismSimulator",rev="for_rmg")); using ReactionMechanismSimulator'
-
- name: Checkout gh-pages Branch
uses: actions/checkout@v2
with:
diff --git a/.gitignore b/.gitignore
index d52d5c680f..38d73e2bc1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -47,7 +47,7 @@ nbproject/*
.vscode/*
# Unit test files
-.coverage
+.coverage*
testing/*
htmlcov/*
@@ -103,6 +103,7 @@ test/rmgpy/test_data/temp_dir_for_testing/cantera/chem001.yaml
rmgpy/test_data/copied_kinetic_lib/
testing/qm/*
test_log.txt
+rmgpy/tools/data/flux/flux/1/*.dot
# example directory - save the inputs but not the outputs
# cantera input files
@@ -117,3 +118,4 @@ examples/**/dictionary.txt
examples/**/reactions.py
examples/**/RMG_libraries
examples/**/species
+examples/**/plots
diff --git a/Dockerfile b/Dockerfile
index cafc18c2a9..5475205334 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -13,24 +13,20 @@ RUN ln -snf /bin/bash /bin/sh
# - libxrender1 required by RDKit
RUN apt-get update && \
apt-get install -y \
- make \
- gcc \
- wget \
- git \
- g++ \
- libxrender1 && \
+ make \
+ gcc \
+ wget \
+ git \
+ g++ \
+ libxrender1 && \
apt-get autoremove -y && \
apt-get clean -y
# Install conda
-RUN wget "https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-Linux-x86_64.sh" && \
- bash Miniforge3-Linux-x86_64.sh -b -p /miniforge && \
- rm Miniforge3-Linux-x86_64.sh
-ENV PATH="$PATH:/miniforge/bin"
-
-# Set solver backend to mamba for speed
-RUN conda install -n base conda-libmamba-solver && \
- conda config --set solver libmamba
+RUN wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh && \
+ bash Miniconda3-latest-Linux-x86_64.sh -b -p /miniconda && \
+ rm Miniconda3-latest-Linux-x86_64.sh
+ENV PATH="$PATH:/miniconda/bin"
# Set Bash as the default shell for following commands
SHELL ["/bin/bash", "-c"]
@@ -50,8 +46,7 @@ RUN git clone --single-branch --branch ${RMG_Py_Branch} --depth 1 https://github
WORKDIR /rmg/RMG-Py
# build the conda environment
-RUN conda env create --file environment.yml && \
- conda clean --all --yes
+RUN conda env create --file environment.yml
# This runs all subsequent commands inside the rmg_env conda environment
#
@@ -60,6 +55,10 @@ RUN conda env create --file environment.yml && \
# in a Dockerfile build script)
SHELL ["conda", "run", "--no-capture-output", "-n", "rmg_env", "/bin/bash", "-c"]
+RUN conda install -c conda-forge julia=1.9.1 pyjulia>=0.6 && \
+ conda install -c rmg pyrms diffeqpy && \
+ conda clean --all --yes
+
# Set environment variables as directed in the RMG installation instructions
ENV RUNNER_CWD=/rmg
ENV PYTHONPATH="$RUNNER_CWD/RMG-Py:$PYTHONPATH"
@@ -70,16 +69,15 @@ ENV PATH="$RUNNER_CWD/RMG-Py:$PATH"
# setting this env variable fixes an issue with Julia precompilation on Windows
ENV JULIA_CPU_TARGET="x86-64,haswell,skylake,broadwell,znver1,znver2,znver3,cascadelake,icelake-client,cooperlake,generic"
RUN make && \
- julia -e 'using Pkg; Pkg.add(PackageSpec(name="PyCall",rev="master")); Pkg.add(PackageSpec(name="ReactionMechanismSimulator",rev=ENV["rmsbranch"])); using ReactionMechanismSimulator' && \
- python -c "import julia; julia.install(); import diffeqpy; diffeqpy.install()"
+ julia -e 'ENV["JULIA_CONDAPKG_BACKEND"] = "Current"; using Pkg; Pkg.add(Pkg.PackageSpec(name="ReactionMechanismSimulator", url="https://github.com/hwpang/ReactionMechanismSimulator.jl.git", rev="fix_installation")); using ReactionMechanismSimulator'
# RMG-Py should now be installed and ready - trigger precompilation and test run
-RUN python-jl rmg.py examples/rmg/minimal/input.py
+RUN python rmg.py examples/rmg/minimal/input.py
# delete the results, preserve input.py
RUN mv examples/rmg/minimal/input.py . && \
rm -rf examples/rmg/minimal/* && \
mv input.py examples/rmg/minimal/
# when running this image, open an interactive bash terminal inside the conda environment
-RUN echo "source activate rmg_env" > ~/.bashrc
+RUN echo "source activate rmg_env" >~/.bashrc
ENTRYPOINT ["/bin/bash", "--login"]
diff --git a/Makefile b/Makefile
index a2f6b3fed2..06055daa3a 100644
--- a/Makefile
+++ b/Makefile
@@ -62,16 +62,16 @@ decython:
find . -name *.pyc -exec rm -f '{}' \;
test-all:
- python-jl -m pytest
+ python -m pytest
test test-unittests:
- python-jl -m pytest -m "not functional and not database"
+ python -m pytest -m "not functional and not database"
test-functional:
- python-jl -m pytest -m "functional"
+ python -m pytest -m "functional"
test-database:
- python-jl -m pytest -m "database"
+ python -m pytest -m "database"
eg0: all
mkdir -p testing/eg0
diff --git a/arkane/encorr/isodesmic.py b/arkane/encorr/isodesmic.py
index 83a9164e17..f4f39d0d85 100644
--- a/arkane/encorr/isodesmic.py
+++ b/arkane/encorr/isodesmic.py
@@ -4,7 +4,7 @@
# #
# RMG - Reaction Mechanism Generator #
# #
-# Copyright (c) 2002-2023 Prof. William H. Green (whgreen@mit.edu), #
+# Copyright (c) 2002-2024 Prof. William H. Green (whgreen@mit.edu), #
# Prof. Richard H. West (r.west@neu.edu) and the RMG Team (rmg_dev@mit.edu) #
# #
# Permission is hereby granted, free of charge, to any person obtaining a #
@@ -41,28 +41,29 @@
https://doi.org/10.1021/jp404158v
"""
-import signal
-from collections import deque
+import logging
+from copy import deepcopy
+from typing import List, Union
-from lpsolve55 import lpsolve, EQ, LE
import numpy as np
-
-from rmgpy.molecule import Molecule
-from rmgpy.quantity import ScalarQuantity
+from scipy.optimize import LinearConstraint, milp
from arkane.modelchem import LOT
-
-# Optional Imports
-try:
- import pyomo.environ as pyo
-except ImportError:
- pyo = None
+from rmgpy.molecule import Bond, Molecule
+from rmgpy.quantity import ScalarQuantity
class ErrorCancelingSpecies:
"""Class for target and known (reference) species participating in an error canceling reaction"""
- def __init__(self, molecule, low_level_hf298, level_of_theory, high_level_hf298=None, source=None):
+ def __init__(
+ self,
+ molecule,
+ low_level_hf298,
+ level_of_theory,
+ high_level_hf298=None,
+ source=None,
+ ):
"""
Args:
@@ -76,35 +77,45 @@ def __init__(self, molecule, low_level_hf298, level_of_theory, high_level_hf298=
if isinstance(molecule, Molecule):
self.molecule = molecule
else:
- raise ValueError(f'ErrorCancelingSpecies molecule attribute must be an rmgpy Molecule object. Instead a '
- f'{type(molecule)} object was given')
+ raise ValueError(
+ f"ErrorCancelingSpecies molecule attribute must be an rmgpy Molecule object. Instead a "
+ f"{type(molecule)} object was given"
+ )
if isinstance(level_of_theory, LOT):
self.level_of_theory = level_of_theory
else:
- raise ValueError(f'The level of theory used to calculate the low level Hf298 must be provided '
- f'for consistency checks. Instead, a {type(level_of_theory)} object was given')
+ raise ValueError(
+ f"The level of theory used to calculate the low level Hf298 must be provided "
+ f"for consistency checks. Instead, a {type(level_of_theory)} object was given"
+ )
if not isinstance(low_level_hf298, ScalarQuantity):
if isinstance(low_level_hf298, tuple):
low_level_hf298 = ScalarQuantity(*low_level_hf298)
else:
- raise TypeError(f'Low level Hf298 should be a ScalarQuantity object or its tuple representation, but '
- f'received {low_level_hf298} instead.')
+ raise TypeError(
+ f"Low level Hf298 should be a ScalarQuantity object or its tuple representation, but "
+ f"received {low_level_hf298} instead."
+ )
self.low_level_hf298 = low_level_hf298
# If the species is a reference species, then the high level data is already known
- if high_level_hf298 is not None and not isinstance(high_level_hf298, ScalarQuantity):
+ if high_level_hf298 is not None and not isinstance(
+ high_level_hf298, ScalarQuantity
+ ):
if isinstance(high_level_hf298, tuple):
high_level_hf298 = ScalarQuantity(*high_level_hf298)
else:
- raise TypeError(f'High level Hf298 should be a ScalarQuantity object or its tuple representation, but '
- f'received {high_level_hf298} instead.')
+ raise TypeError(
+ f"High level Hf298 should be a ScalarQuantity object or its tuple representation, but "
+ f"received {high_level_hf298} instead."
+ )
self.high_level_hf298 = high_level_hf298
self.source = source
def __repr__(self):
- return f''
+ return f""
class ErrorCancelingReaction:
@@ -131,22 +142,24 @@ def __init__(self, target, species):
# Perform a consistency check that all species are using the same level of theory
for spcs in species.keys():
if spcs.level_of_theory != self.level_of_theory:
- raise ValueError(f'Species {spcs} has level of theory {spcs.level_of_theory}, which does not match the '
- f'level of theory of the reaction of {self.level_of_theory}')
+ raise ValueError(
+ f"Species {spcs} has level of theory {spcs.level_of_theory}, which does not match the "
+ f"level of theory of the reaction of {self.level_of_theory}"
+ )
# Does not include the target, which is handled separately.
self.species = species
def __repr__(self):
- reactant_string = f'1*{self.target.molecule.to_smiles()}'
- product_string = ''
+ reactant_string = f"1*{self.target.molecule.to_smiles()}"
+ product_string = ""
for spcs, coeff in self.species.items():
if coeff > 0:
- product_string += f' + {int(coeff)}*{spcs.molecule.to_smiles()}'
+ product_string += f" + {int(coeff)}*{spcs.molecule.to_smiles()}"
else:
- reactant_string += f' + {-1*int(coeff)}*{spcs.molecule.to_smiles()}'
+ reactant_string += f" + {-1*int(coeff)}*{spcs.molecule.to_smiles()}"
- return f' {product_string[3:]} >'
+ return f" {product_string[3:]} >"
def calculate_target_thermo(self):
"""
@@ -155,12 +168,191 @@ def calculate_target_thermo(self):
Returns:
rmgpy.quantity.ScalarQuantity: Hf298 in 'J/mol' estimated for the target species
"""
- low_level_h_rxn = sum(spec[0].low_level_hf298.value_si*spec[1] for spec in self.species.items()) - \
- self.target.low_level_hf298.value_si
+ low_level_h_rxn = (
+ sum(
+ spec.low_level_hf298.value_si * coeff
+ for spec, coeff in self.species.items()
+ )
+ - self.target.low_level_hf298.value_si
+ )
+
+ target_thermo = (
+ sum(
+ spec.high_level_hf298.value_si * coeff
+ for spec, coeff in self.species.items()
+ )
+ - low_level_h_rxn
+ )
+ return ScalarQuantity(target_thermo, "J/mol")
+
+
+class AtomConstraint:
+ def __init__(self, label, connections=None):
+ self.label = label
+ self.connections = connections if connections is not None else []
+
+ def __eq__(self, other):
+ if isinstance(other, AtomConstraint):
+ if self.label == other.label:
+ return self.match_connections(other)
+
+ return False
+
+ else:
+ raise NotImplementedError(
+ f"AtomConstraint object has no __eq__ defined for other object of type "
+ f"{type(other)}"
+ )
+
+ def match_connections(self, other):
+ if len(self.connections) != len(other.connections):
+ return False
+
+ connections = deepcopy(other.connections)
+ for c in self.connections:
+ for i, c_other in enumerate(connections):
+ if c == c_other:
+ break
+ else:
+ return False
+ connections.pop(i)
+
+ return True
+
+ def __repr__(self):
+ return f"{self.label}" + "".join([f"({c})" for c in self.connections])
+
+
+class BondConstraint:
+ def __init__(self, atom1, atom2, bond_order):
+ self.atom1 = atom1
+ self.atom2 = atom2
+ self.bond_order = bond_order
+
+ def __eq__(self, other):
+ if isinstance(other, BondConstraint):
+ if self.bond_order == other.bond_order:
+ if (self.atom1 == other.atom1 and self.atom2 == other.atom2) or (
+ self.atom1 == other.atom2 and self.atom2 == other.atom1
+ ):
+ return True
+ return False
+
+ if isinstance(other, GenericConstraint):
+ return False
+
+ else:
+ raise NotImplementedError(
+ f"BondConstraint object has no __eq__ defined for other object of type "
+ f"{type(other)}"
+ )
+
+ def __repr__(self):
+ symbols = ["", "-", "=", "#"]
+ return f"{self.atom1}{symbols[self.bond_order]}{self.atom2}"
+
+
+class Connection:
+ def __init__(self, atom, bond_order):
+ self.atom = atom
+ self.bond_order = bond_order
+
+ def __eq__(self, other):
+ if isinstance(other, Connection):
+ if self.bond_order == other.bond_order:
+ if self.atom == other.atom:
+ return True
+ return False
+
+ else:
+ raise NotImplementedError(
+ f"Connection object has no __eq__ defined for other object of type {type(other)}"
+ )
+
+ def __repr__(self):
+ symbols = ["", "-", "=", "#"]
+ return f"{symbols[self.bond_order]}{self.atom}"
+
+
+class GenericConstraint:
+ def __init__(self, constraint_str: str):
+ self.constraint_str = constraint_str
+
+ def __eq__(self, other):
+ if isinstance(other, BondConstraint):
+ return False
+ elif isinstance(other, GenericConstraint):
+ return self.constraint_str == other.constraint_str
+ else:
+ raise NotImplementedError(
+ f"GenericConstraint object has no __eq__ defined for other object of "
+ f"type {type(other)}"
+ )
+ def __repr__(self):
+ return self.constraint_str
+
+
+def bond_centric_constraints(
+ species: ErrorCancelingSpecies, constraint_class: str
+) -> List[BondConstraint]:
+ constraints = []
+ contraint_func = CONSTRAINT_CLASSES[constraint_class]
+ molecule = species.molecule
+
+ for bond in molecule.get_all_edges():
+ constraints.append(contraint_func(bond))
+
+ return constraints
+
+
+def _buerger_rc2(bond: Bond) -> BondConstraint:
+ atom1 = bond.atom1
+ atom2 = bond.atom2
+
+ if atom1.symbol > atom2.symbol:
+ atom1, atom2 = atom2, atom1
+
+ atom1 = AtomConstraint(label=atom1.symbol)
+ atom2 = AtomConstraint(label=atom2.symbol)
+
+ return BondConstraint(atom1=atom1, atom2=atom2, bond_order=int(bond.order))
+
+
+def _buerger_rc3(bond: Bond) -> BondConstraint:
+ atom1 = bond.atom1
+ atom2 = bond.atom2
+
+ if atom1.symbol > atom2.symbol:
+ atom1, atom2 = atom2, atom1
- target_thermo = sum(spec[0].high_level_hf298.value_si*spec[1] for spec in self.species.items()) - \
- low_level_h_rxn
- return ScalarQuantity(target_thermo, 'J/mol')
+ atom1 = AtomConstraint(label=f"{atom1.symbol}{atom1.connectivity1}")
+ atom2 = AtomConstraint(label=f"{atom2.symbol}{atom2.connectivity1}")
+
+ return BondConstraint(atom1=atom1, atom2=atom2, bond_order=int(bond.order))
+
+
+def _buerger_rc4(bond: Bond) -> BondConstraint:
+ atom1 = bond.atom1
+ atom2 = bond.atom2
+
+ if atom1.symbol > atom2.symbol:
+ atom1, atom2 = atom2, atom1
+
+ atoms = []
+
+ for atom in [atom1, atom2]:
+ connections = []
+ for a, b in atom.bonds.items():
+ ac = AtomConstraint(label=f"{a.symbol}{a.connectivity1}")
+ bond_order = b.order
+ connections.append(Connection(atom=ac, bond_order=bond_order))
+ atoms.append(
+ AtomConstraint(
+ label=f"{atom.symbol}{atom.connectivity1}", connections=connections
+ )
+ )
+
+ return BondConstraint(atom1=atoms[0], atom2=atoms[1], bond_order=int(bond.order))
class SpeciesConstraints:
@@ -168,76 +360,181 @@ class SpeciesConstraints:
A class for defining and enumerating constraints to ReferenceSpecies objects for error canceling reactions
"""
- def __init__(self, target, reference_list, conserve_bonds=True, conserve_ring_size=True):
+ def __init__(
+ self,
+ target,
+ reference_list,
+ isodesmic_class="rc2",
+ conserve_ring_size=True,
+ limit_charges=True,
+ limit_scope=True,
+ ):
"""
Define the constraints that will be enforced, and determine the mapping of indices in the constraint vector to
the labels for these constraints.
- To reduce the size of the linear programming problem that will try to find error canceling reactions of the
- target and subsets of the reference species, the `reference_species` list is automatically pruned to remove
- species that have additional atom, bond, and/or ring attributes not found in the target molecule.
+ Notes:
+ To reduce the size of the linear programming problem that will try to find error canceling reactions of the
+ target and subsets of the reference species, the `reference_species` list is automatically pruned to remove
+ species that have additional atom, bond, and/or ring attributes not found in the target molecule.
+
+ Charge is also explicitly conserved, as there are charged species in the reference database
Args:
target (ErrorCancelingSpecies): The target species of the error canceling reaction scheme
reference_list(list): A list of ErrorCancelingSpecies objects for the reference
species that can participate in the error canceling reaction scheme
- conserve_bonds (bool, optional): Enforce the number of each bond type be conserved
- conserve_ring_size (bool, optional): Enforce that the number of each ring size be conserved
+ isodesmic_class (str, optional): Reaction classes as defined by Buerger et al. that determine how specific
+ the constraints are.
+ conserve_ring_size (bool, optional): Enforce that the number of each ring size be conserved.
+ limit_charges (bool, optional): Only allow species in the reaction that are within the range [C, 0] for
+ anions or [0, C] for cations where "C" is the charge of the target
+ limit_scope (bool, optional): Exclude any molecules from the reference set that have features not contained
+ in the target molecule. This will reduce the size of the linear programing problem being solved to yield
+ faster, possibly more accurate results
"""
self.target = target
self.all_reference_species = reference_list
self.reference_species = []
- self.conserve_bonds = conserve_bonds
+ self.isodesmic_class = isodesmic_class
self.conserve_ring_size = conserve_ring_size
- self.constraint_map = self._get_constraint_map()
-
- def _get_constraint_map(self):
- # Enumerate all of the constraints in the target molecule to initialize the constraint mapping
- constraint_map = {label: i for i, label in enumerate(self.target.molecule.get_element_count().keys())}
- if self.conserve_bonds:
- j = len(constraint_map)
- constraint_map.update(
- {label: j + i for i, label in enumerate(self.target.molecule.enumerate_bonds().keys())})
+ self.limit_charges = limit_charges
+ self.limit_scope = limit_scope
+
+ def _get_constraint_lists(self):
+ full_constraints_list = [self._get_all_constraints(self.target)]
+ for ref_spcs in self.all_reference_species:
+ full_constraints_list.append(self._get_all_constraints(ref_spcs))
+
+ return full_constraints_list
+
+ def _get_ring_constraints(
+ self, species: ErrorCancelingSpecies
+ ) -> List[GenericConstraint]:
+ ring_features = []
+ rings = species.molecule.get_smallest_set_of_smallest_rings()
+ for ring in rings:
+ ring_features.append(GenericConstraint(constraint_str=f"{len(ring)}_ring"))
+
+ return ring_features
+
+ def _get_all_constraints(
+ self, species: ErrorCancelingSpecies
+ ) -> List[Union[BondConstraint, GenericConstraint]]:
+ features = bond_centric_constraints(species, self.isodesmic_class)
if self.conserve_ring_size:
- j = len(constraint_map)
- possible_rings_sizes = set(f'{len(sssr)}_ring' for sssr in
- self.target.molecule.get_smallest_set_of_smallest_rings())
- constraint_map.update({label: j + i for i, label in enumerate(possible_rings_sizes)})
+ features += self._get_ring_constraints(species)
- return constraint_map
+ features.sort(key=lambda x: x.__repr__())
+ return features
- def _enumerate_constraints(self, species):
+ def _enumerate_constraints(self, full_constraints_list):
+ """
+ Return the target constraint counts and the reference constraint counts.
"""
- Determine the constraint vector for a molecule given the enforced constraints
- Args:
- species (ErrorCancelingSpecies): Species whose constraints are to be enumerated
+ target_constraints = full_constraints_list[0]
+ reference_constraintss = full_constraints_list[1:]
- Returns:
- np.ndarray: vector of the number of instances of each constraining feature e.g. number of carbon atoms
- """
- constraint_vector = np.zeros(len(self.constraint_map))
- molecule = species.molecule
+ # Enumerate through the constraints of reference species and keep only those that are present in the target
+ enumerated_reference_constraintss = []
+
+ if self.limit_scope:
+
+ for i, ref_spc in enumerate(self.all_reference_species):
- try:
- atoms = molecule.get_element_count()
- for atom_label, count in atoms.items():
- constraint_vector[self.constraint_map[atom_label]] += count
+ if not all(constraint in target_constraints for constraint in reference_constraintss[i]):
+ continue
- if self.conserve_bonds:
- bonds = molecule.enumerate_bonds()
- for bond_label, count in bonds.items():
- constraint_vector[self.constraint_map[bond_label]] += count
+ self.reference_species.append(ref_spc)
+ enumerated_reference_constraintss.append(reference_constraintss[i])
- if self.conserve_ring_size:
- rings = molecule.get_smallest_set_of_smallest_rings()
- for ring in rings:
- constraint_vector[self.constraint_map[f'{len(ring)}_ring']] += 1
- except KeyError: # This molecule has a feature not found in the target molecule. Return None to exclude this
- return None
+ else:
+ self.reference_species = self.all_reference_species
+ enumerated_reference_constraintss = reference_constraintss
+
+ # Get a list of the unique constraints sorted by their string representation
+ if self.limit_scope:
- return constraint_vector
+ # The scope of constraints to consider is the target constraints
+ unique_constraints = self._get_unique_constraints(target_constraints)
+ unique_constraints.sort(key=lambda x: x.__repr__())
+
+ else:
+ all_constraints = target_constraints + [constraint for constraints in enumerated_reference_constraintss for constraint in constraints]
+ unique_constraints = self._get_unique_constraints(all_constraints)
+ unique_constraints.sort(key=lambda x: x.__repr__())
+
+ # Get the counts of each unique constraint in the target and reference constraints
+ target_constraint_counts = [target_constraints.count(c) for c in unique_constraints]
+ reference_constraint_counts = []
+
+ for i, ref_constraints in enumerate(enumerated_reference_constraintss):
+ reference_constraint_counts.append([ref_constraints.count(c) for c in unique_constraints])
+
+ return target_constraint_counts, reference_constraint_counts
+
+ def _get_unique_constraints(self, constraints):
+ # Constraints are unhashable, so we need to use some workaround to get unique constraints
+ constraint_dict = {constraint.__repr__(): constraint for constraint in constraints}
+ return list(constraint_dict.values())
+
+ def _enumerate_charge_constraints(self, target_constraints, reference_constraints):
+ charge = self.target.molecule.get_net_charge()
+ target_constraints.append(charge)
+
+ for i, spcs in enumerate(self.reference_species):
+ reference_constraints[i].append(spcs.molecule.get_net_charge())
+
+ if self.limit_charges:
+ allowed_reference_species = []
+ new_constraints = []
+
+ if charge < 0:
+ allowable_charges = list(range(charge, 0))
+ elif charge > 0:
+ allowable_charges = list(range(1, charge + 1))
+ else:
+ allowable_charges = [0]
+ for i, spcs in enumerate(self.reference_species):
+ if reference_constraints[i][-1] in allowable_charges:
+ allowed_reference_species.append(spcs)
+ new_constraints.append(reference_constraints[i])
+
+ reference_constraints = new_constraints
+ self.reference_species = allowed_reference_species
+
+ return target_constraints, reference_constraints
+
+ def _enumerate_element_constraints(self, target_constraints, reference_constraints):
+ all_elements = set()
+ for spc in self.reference_species:
+ all_elements.update(spc.molecule.get_element_count().keys())
+
+ # Check that the target and reference species have the same elements to be able to satisfy mass conservation
+ if set(self.target.molecule.get_element_count().keys()) != all_elements:
+ logging.warning(
+ "Target species and reference species do not have the same elements:"
+ f"Target: {' '.join(self.target.molecule.get_element_count().keys())}"
+ f"Reference: {all_elements}"
+ )
+
+ all_elements.update(self.target.molecule.get_element_count().keys())
+ all_elements = sorted(list(all_elements))
+
+ element_count = self.target.molecule.get_element_count()
+ new_constraints = [element_count.get(element, 0) for element in all_elements]
+ target_constraints.extend(new_constraints)
+
+ for i, spc in enumerate(self.reference_species):
+ element_count = spc.molecule.get_element_count()
+ new_constraints = [
+ element_count.get(element, 0) for element in all_elements
+ ]
+ reference_constraints[i].extend(new_constraints)
+
+ return target_constraints, reference_constraints
def calculate_constraints(self):
"""
@@ -248,16 +545,41 @@ def calculate_constraints(self):
- target constraint vector (1 x len(constraints))
- constraint matrix for allowable reference species (len(self.reference_species) x len(constraints))
"""
- target_constraints = self._enumerate_constraints(self.target)
- constraint_matrix = []
-
- for spcs in self.all_reference_species:
- spcs_constraints = self._enumerate_constraints(spcs)
- if spcs_constraints is not None: # This species is allowed
- self.reference_species.append(spcs)
- constraint_matrix.append(spcs_constraints)
+ full_constraint_list = self._get_constraint_lists()
+ target_constraints, reference_constraints = self._enumerate_constraints(
+ full_constraint_list
+ )
+ target_constraints, reference_constraints = self._enumerate_charge_constraints(
+ target_constraints, reference_constraints
+ )
+ target_constraints, reference_constraints = self._enumerate_element_constraints(
+ target_constraints, reference_constraints
+ )
+
+ target_constraints = np.array(target_constraints, dtype=int)
+ constraint_matrix = np.array(reference_constraints, dtype=int)
+
+ return target_constraints, constraint_matrix
+
+
+def _clean_up_constraints(target_constraints, constraint_matrix):
+ # make sure that the constraint matrix is 2d
+ if len(constraint_matrix.shape) == 1:
+ constraint_matrix = np.array([constraint_matrix], dtype=int)
+
+ # Remove any columns that are all zeros
+ zero_indices = np.where(~constraint_matrix.any(axis=0))[0]
+ # Check that this wouldn't eliminate a non-zero target entry
+ for z in zero_indices:
+ if (
+ target_constraints[z] != 0
+ ): # This problem is not solvable. Return None to signal this
+ return None, None
+ indices = [i for i in range(constraint_matrix.shape[1]) if i not in zero_indices]
+ constraint_matrix = np.take(constraint_matrix, indices=indices, axis=1)
+ target_constraints = np.take(target_constraints, indices=indices)
- return target_constraints, np.array(constraint_matrix, dtype=int)
+ return target_constraints, constraint_matrix
class ErrorCancelingScheme:
@@ -265,27 +587,42 @@ class ErrorCancelingScheme:
A Base class for calculating target species thermochemistry using error canceling reactions
"""
- def __init__(self, target, reference_set, conserve_bonds, conserve_ring_size):
+ def __init__(
+ self,
+ target,
+ reference_set,
+ isodesmic_class,
+ conserve_ring_size,
+ limit_charges,
+ limit_scope,
+ ):
"""
Args:
target (ErrorCancelingSpecies): Species whose Hf298 will be calculated using error canceling reactions
reference_set (list): list of reference species (as ErrorCancelingSpecies) that can participate
in error canceling reactions to help calculate the thermochemistry of the target
- conserve_bonds (bool): Flag to determine if the number and type of each bond must be conserved in each error
- canceling reaction
conserve_ring_size (bool): Flag to determine if the number of each ring size must be conserved in each error
canceling reaction
"""
self.target = target
- self.constraints = SpeciesConstraints(target, reference_set, conserve_bonds=conserve_bonds,
- conserve_ring_size=conserve_ring_size)
-
- self.target_constraint, self.constraint_matrix = self.constraints.calculate_constraints()
+ self.constraints = SpeciesConstraints(
+ target,
+ reference_set,
+ isodesmic_class=isodesmic_class,
+ conserve_ring_size=conserve_ring_size,
+ limit_charges=limit_charges,
+ limit_scope=limit_scope,
+ )
+
+ (
+ self.target_constraint,
+ self.constraint_matrix,
+ ) = self.constraints.calculate_constraints()
self.reference_species = self.constraints.reference_species
- def _find_error_canceling_reaction(self, reference_subset, milp_software=None):
+ def _find_error_canceling_reaction(self, reference_subset):
"""
Automatically find a valid error canceling reaction given a subset of the available benchmark species. This
is done by solving a mixed integer linear programming (MILP) problem similiar to
@@ -293,107 +630,45 @@ def _find_error_canceling_reaction(self, reference_subset, milp_software=None):
Args:
reference_subset (list): A list of indices from self.reference_species that can participate in the reaction
- milp_software (list, optional): Solvers to try in order. Defaults to ['lpsolve'] or if pyomo is available
- defaults to ['lpsolve', 'pyomo']. lpsolve is usually faster.
Returns:
tuple(ErrorCancelingReaction, np.ndarray)
- Reaction with the target species (if a valid reaction is found, else ``None``)
- Indices (of the subset) for the species that participated in the return reaction
"""
- if milp_software is None:
- milp_software = ['lpsolve']
- if pyo is not None:
- milp_software.append('pyomo')
# Define the constraints based on the provided subset
c_matrix = np.take(self.constraint_matrix, reference_subset, axis=0)
+
+ # Remove unnecessary constraints
+ target_constraint, c_matrix = _clean_up_constraints(
+ self.target_constraint, c_matrix
+ )
+ if target_constraint is None or c_matrix is None: # The problem is not solvable
+ return None, None
+
+ # Setup MILP problem
c_matrix = np.tile(c_matrix, (2, 1))
sum_constraints = np.sum(c_matrix, 1, dtype=int)
- targets = -1*self.target_constraint
+ targets = -1 * target_constraint
m = c_matrix.shape[0]
n = c_matrix.shape[1]
- split = int(m/2)
-
- for solver in milp_software:
- if solver == 'pyomo':
- # Check that pyomo is available
- if pyo is None:
- raise ImportError('Cannot import optional package pyomo. Either install this dependency with '
- '`conda install -c conda-forge pyomo glpk` or set milp_software to `lpsolve`')
-
- # Setup the MILP problem using pyomo
- lp_model = pyo.ConcreteModel()
- lp_model.i = pyo.RangeSet(0, m - 1)
- lp_model.j = pyo.RangeSet(0, n - 1)
- lp_model.r = pyo.RangeSet(0, split-1) # indices before the split correspond to reactants
- lp_model.p = pyo.RangeSet(split, m - 1) # indices after the split correspond to products
- lp_model.v = pyo.Var(lp_model.i, domain=pyo.NonNegativeIntegers) # The stoich. coef. we are solving for
- lp_model.c = pyo.Param(lp_model.i, lp_model.j, initialize=lambda _, i_ind, j_ind: c_matrix[i_ind,
- j_ind])
- lp_model.s = pyo.Param(lp_model.i, initialize=lambda _, i_ind: sum_constraints[i_ind])
- lp_model.t = pyo.Param(lp_model.j, initialize=lambda _, j_ind: targets[j_ind])
-
- lp_model.obj = pyo.Objective(rule=_pyo_obj_expression)
- lp_model.constraints = pyo.Constraint(lp_model.j, rule=_pyo_constraint_rule)
-
- # Solve the MILP problem using the GLPK MILP solver (https://www.gnu.org/software/glpk/)
- opt = pyo.SolverFactory('glpk')
- results = opt.solve(lp_model, timelimit=1)
-
- # Return the solution if a valid reaction is found. Otherwise continue to next solver
- if results.solver.termination_condition == pyo.TerminationCondition.optimal:
- # Extract the solution and find the species with non-zero stoichiometric coefficients
- solution = lp_model.v.extract_values().values()
- break
+ split = int(m / 2)
- elif solver == 'lpsolve':
- # Save the current signal handler
- sig = signal.getsignal(signal.SIGINT)
-
- # Setup the MILP problem using lpsolve
- lp = lpsolve('make_lp', 0, m)
- lpsolve('set_verbose', lp, 2) # Reduce the logging from lpsolve
- lpsolve('set_obj_fn', lp, sum_constraints)
- lpsolve('set_minim', lp)
-
- for j in range(n):
- lpsolve('add_constraint', lp, np.concatenate((c_matrix[:split, j], -1*c_matrix[split:, j])), EQ,
- targets[j])
-
- lpsolve('add_constraint', lp, np.ones(m), LE, 20) # Use at most 20 species (including replicates)
- lpsolve('set_timeout', lp, 1) # Move on if lpsolve can't find a solution quickly
-
- # Constrain v_i to be 4 or less
- for i in range(m):
- lpsolve('set_upbo', lp, i, 4)
-
- # All v_i must be integers
- lpsolve('set_int', lp, [True]*m)
-
- status = lpsolve('solve', lp)
-
- # Reset signal handling since lpsolve changed it
- try:
- signal.signal(signal.SIGINT, sig)
- except ValueError:
- # This is not being run in the main thread, so we cannot reset signal
- pass
- except TypeError:
- print(
- "Failed to reset signal handling in LPSolve - are you running pytest?"
- )
+ constraints = tuple((LinearConstraint(A=np.concatenate((c_matrix[:split, j], -1 * c_matrix[split:, j])), lb=targets[j], ub=targets[j]) for j in range(n)))
- # Return the solution if a valid reaction is found. Otherwise continue to next solver
- if status == 0:
- _, solution = lpsolve('get_solution', lp)[:2]
- break
+ result = milp(
+ sum_constraints,
+ integrality=1,
+ constraints=constraints,
+ options={"time_limit": 10},
+ )
- else:
- raise ValueError(f'Unrecognized MILP solver {solver} for isodesmic reaction generation')
-
- else:
+ if result.status != 0:
+ logging.warning("Optimization could not find a valid solution.")
return None, None
+
+ solution = result.x
reaction = ErrorCancelingReaction(self.target, dict())
subset_indices = []
@@ -401,13 +676,19 @@ def _find_error_canceling_reaction(self, reference_subset, milp_software=None):
if v > 0:
subset_indices.append(index % split)
if index < split:
- reaction.species.update({self.reference_species[reference_subset[index]]: -v})
+ reaction.species.update(
+ {self.reference_species[reference_subset[index]]: -v}
+ )
else:
- reaction.species.update({self.reference_species[reference_subset[index % split]]: v})
+ reaction.species.update(
+ {self.reference_species[reference_subset[index % split]]: v}
+ )
return reaction, np.array(subset_indices)
- def multiple_error_canceling_reaction_search(self, n_reactions_max=20, milp_software=None):
+ def multiple_error_canceling_reaction_search(
+ self, n_reactions_max=20,
+ ):
"""
Generate multiple error canceling reactions involving the target and a subset of the reference species.
@@ -418,21 +699,20 @@ def multiple_error_canceling_reaction_search(self, n_reactions_max=20, milp_soft
Args:
n_reactions_max (int, optional): The maximum number of found reactions that will be returned, after which no
further searching will occur even if there are possible subsets left in the queue.
- milp_software (list, optional): Solvers to try in order. Defaults to ['lpsolve'] or if pyomo is available
- defaults to ['lpsolve', 'pyomo']. lpsolve is usually faster.
Returns:
list: A list of the found error canceling reactions
"""
- subset_queue = deque()
- subset_queue.append(np.arange(0, len(self.reference_species)))
+ subset_queue = [np.arange(0, len(self.reference_species))]
reaction_list = []
while (len(subset_queue) != 0) and (len(reaction_list) < n_reactions_max):
- subset = subset_queue.popleft()
+ subset = subset_queue.pop()
if len(subset) == 0:
continue
- reaction, subset_indices = self._find_error_canceling_reaction(subset, milp_software=milp_software)
+ reaction, subset_indices = self._find_error_canceling_reaction(
+ subset
+ )
if reaction is None:
continue
else:
@@ -441,9 +721,18 @@ def multiple_error_canceling_reaction_search(self, n_reactions_max=20, milp_soft
for index in subset_indices:
subset_queue.append(np.delete(subset, index))
+ # Clean up the queue to remove subsets that would allow for already found reactions
+ new_queue = []
+ reaction_indices = [subset[i] for i in subset_indices]
+ for s in subset_queue:
+ if not all([i in s for i in reaction_indices]):
+ new_queue.append(s)
+
+ subset_queue = new_queue
+
return reaction_list
- def calculate_target_enthalpy(self, n_reactions_max=20, milp_software=None):
+ def calculate_target_enthalpy(self, n_reactions_max=5):
"""
Perform a multiple error canceling reactions search and calculate hf298 for the target species by taking the
median hf298 value from among the error canceling reactions found
@@ -451,47 +740,57 @@ def calculate_target_enthalpy(self, n_reactions_max=20, milp_software=None):
Args:
n_reactions_max (int, optional): The maximum number of found reactions that will returned, after which no
further searching will occur even if there are possible subsets left in the queue.
- milp_software (list, optional): Solvers to try in order. Defaults to ['lpsolve'] or if pyomo is available
- defaults to ['lpsolve', 'pyomo']. lpsolve is usually faster.
Returns:
tuple(ScalarQuantity, list)
- Standard heat of formation at 298 K calculated for the target species
- reaction list containing all error canceling reactions found
"""
- reaction_list = self.multiple_error_canceling_reaction_search(n_reactions_max, milp_software)
+ reaction_list = self.multiple_error_canceling_reaction_search(n_reactions_max)
+ if len(reaction_list) == 0: # No reactions found
+ return None, reaction_list
h298_list = np.zeros(len(reaction_list))
for i, rxn in enumerate(reaction_list):
h298_list[i] = rxn.calculate_target_thermo().value_si
- return ScalarQuantity(np.median(h298_list), 'J/mol'), reaction_list
-
-
-def _pyo_obj_expression(model):
- return pyo.summation(model.v, model.s, index=model.i)
-
-
-def _pyo_constraint_rule(model, col):
- return sum(model.v[row] * model.c[row, col] for row in model.r) - \
- sum(model.v[row] * model.c[row, col] for row in model.p) == model.t[col]
+ return ScalarQuantity(np.median(h298_list), "J/mol"), reaction_list
class IsodesmicScheme(ErrorCancelingScheme):
"""
An error canceling reaction where the number and type of both atoms and bonds are conserved
"""
+
def __init__(self, target, reference_set):
- super().__init__(target, reference_set, conserve_bonds=True, conserve_ring_size=False)
+ super().__init__(
+ target,
+ reference_set,
+ isodesmic_class="rc2",
+ conserve_ring_size=False,
+ limit_charges=True,
+ limit_scope=True,
+ )
class IsodesmicRingScheme(ErrorCancelingScheme):
"""
A stricter form of the traditional isodesmic reaction scheme where the number of each ring size is also conserved
"""
+
def __init__(self, target, reference_set):
- super().__init__(target, reference_set, conserve_bonds=True, conserve_ring_size=True)
+ super().__init__(
+ target,
+ reference_set,
+ isodesmic_class="rc2",
+ conserve_ring_size=True,
+ limit_charges=True,
+ limit_scope=True,
+ )
+
+
+CONSTRAINT_CLASSES = {"rc2": _buerger_rc2, "rc3": _buerger_rc3, "rc4": _buerger_rc4}
-if __name__ == '__main__':
+if __name__ == "__main__":
pass
diff --git a/arkane/statmech.py b/arkane/statmech.py
index 4380bad121..fdfac0b512 100644
--- a/arkane/statmech.py
+++ b/arkane/statmech.py
@@ -37,6 +37,7 @@
import math
import os
import pathlib
+import warnings
import matplotlib.pyplot as plt
import numpy as np
@@ -1132,8 +1133,8 @@ def project_rotors(conformer, hessian, rotors, linear, is_ts, get_projected_out_
logging.debug('Frequencies from internal Hessian')
for i in range(3 * n_atoms - external):
- with np.warnings.catch_warnings():
- np.warnings.filterwarnings('ignore', r'invalid value encountered in sqrt')
+ with warnings.catch_warnings():
+ warnings.filterwarnings('ignore', r'invalid value encountered in sqrt')
logging.debug(np.sqrt(eig[i]) / (2 * math.pi * constants.c * 100))
# Now we can start thinking about projecting out the internal rotations
@@ -1243,8 +1244,8 @@ def project_rotors(conformer, hessian, rotors, linear, is_ts, get_projected_out_
logging.debug('Frequencies from projected Hessian')
for i in range(3 * n_atoms):
- with np.warnings.catch_warnings():
- np.warnings.filterwarnings('ignore', r'invalid value encountered in sqrt')
+ with warnings.catch_warnings():
+ warnings.filterwarnings('ignore', r'invalid value encountered in sqrt')
logging.debug(np.sqrt(eig[i]) / (2 * math.pi * constants.c * 100))
return np.sqrt(eig[-n_vib:]) / (2 * math.pi * constants.c * 100)
diff --git a/documentation/Makefile b/documentation/Makefile
index d42279f698..5c30e66cd7 100644
--- a/documentation/Makefile
+++ b/documentation/Makefile
@@ -3,7 +3,7 @@
# You can set these variables from the command line.
SPHINXOPTS =
-SPHINXBUILD = python-jl $$(which sphinx-build)
+SPHINXBUILD = python $$(which sphinx-build)
PAPER =
BUILDDIR = build
diff --git a/documentation/source/users/rmg/installation/anacondaDeveloper.rst b/documentation/source/users/rmg/installation/anacondaDeveloper.rst
index e1aed5f7db..bcfc906039 100644
--- a/documentation/source/users/rmg/installation/anacondaDeveloper.rst
+++ b/documentation/source/users/rmg/installation/anacondaDeveloper.rst
@@ -6,6 +6,11 @@ Installation by Source Using Anaconda Environment for Unix-based Systems: Linux
#. Install the `conda` package manager via `miniforge`, if you do not already have it (or Anaconda), by following the `Miniforge installation instructions `_.
+#. If your `conda` version is older than 23.10.0, switch the solver backend to `libmamba` ::
+
+ conda install -n base conda-libmamba-solver
+ conda config --set solver libmamba
+
#. There are a few system-level dependencies which are required and should not be installed via Conda. These include
`Git `_ for version control, `GNU Make `_, and the C and C++ compilers from the `GNU Compiler Collection (GCC) `_ for compiling RMG.
@@ -52,11 +57,6 @@ Installation by Source Using Anaconda Environment for Unix-based Systems: Linux
For information on using ``ssh`` with GitHub see the `Connecting to GitHub with SSH `_
-#. Switch the conda solver backend to speed up creation of the RMG environment ::
-
- conda install -n base conda-libmamba-solver
- conda config --set solver libmamba
-
#. Navigate to the RMG-Py directory ::
cd RMG-Py
@@ -91,16 +91,11 @@ Installation by Source Using Anaconda Environment for Unix-based Systems: Linux
conda activate rmg_env
-#. Switch the conda solver to libmamba again, to accelerate any changes you might make to this conda environment in the future::
-
- conda config --set solver libmamba
-
#. Compile RMG-Py after activating the conda environment ::
make
#. Modify environment variables. Add RMG-Py to the PYTHONPATH to ensure that you can access RMG modules from any folder.
- *This is important before the next step in which julia dependencies are installed.*
Also, add your RMG-Py folder to PATH to launch ``rmg.py`` from any folder.
In general, these commands should be placed in the appropriate shell initialization file.
@@ -115,16 +110,17 @@ Installation by Source Using Anaconda Environment for Unix-based Systems: Linux
Be sure to either close and reopen your terminal to refresh your environment variables (``source ~/.bashrc`` or ``source ~/.zshrc``).
-#. Install and Link Julia dependencies: ::
+#. **Optional (Recommended)**: Install and Link Julia dependencies. Ensure that you have modified your environment variables as described above, and then run the following: ::
julia -e 'using Pkg; Pkg.add("PyCall");Pkg.build("PyCall");Pkg.add(PackageSpec(name="ReactionMechanismSimulator",rev="for_rmg")); using ReactionMechanismSimulator;'
python -c "import julia; julia.install(); import diffeqpy; diffeqpy.install()"
+ Installing these dependencies will allow using ``method='ode'`` when solving the Master Equation with Arkane and using ``ReactionMechanismSimulator.jl``-based reactors in RMG.
#. Finally, you can run RMG from any location by typing the following (given that you have prepared the input file as ``input.py`` in the current folder). ::
- python-jl replace/with/path/to/rmg.py input.py
+ python replace/with/path/to/rmg.py input.py
You may now use RMG-Py, Arkane, as well as any of the :ref:`Standalone Modules ` included in the RMG-Py package.
@@ -139,7 +135,7 @@ You might have to edit them slightly to match your exact paths. Specifically,
you will need ``/opt/miniconda3/envs/rmg_env`` to point to where your conda environment is located.
This configuration will allow you to debug the rms_constant_V example, running through
-python-jl. ::
+python. ::
{
"name": "Python: rmg.py rms_constant_V",
@@ -147,7 +143,7 @@ python-jl. ::
"request": "launch",
"cwd": "${workspaceFolder}/",
"program": "rmg.py",
- "python": "/opt/miniconda3/envs/rmg_env/bin/python-jl",
+ "python": "/opt/miniconda3/envs/rmg_env/bin/python",
"args": [
"examples/rmg/rms_constant_V/input.py",
],
@@ -166,7 +162,7 @@ Open one of the many test files named ``*Test.py`` in ``test/rmgpy`` before you
"type": "python",
"request": "launch",
"program": "/opt/miniconda3/envs/rmg_env/bin/pytest",
- "python": "/opt/miniconda3/envs/rmg_env/bin/python-jl",
+ "python": "/opt/miniconda3/envs/rmg_env/bin/python",
"args": [
"--capture=no",
"--verbose",
@@ -186,7 +182,7 @@ This configuration will allow you to debug running all the database tests.::
"type": "python",
"request": "launch",
"program": "/opt/miniconda3/envs/rmg_env/bin/pytest",
- "python": "/opt/miniconda3/envs/rmg_env/bin/python-jl",
+ "python": "/opt/miniconda3/envs/rmg_env/bin/python",
"args": [
"--capture=no",
"--verbose",
@@ -207,7 +203,7 @@ This configuration will allow you to use the debugger breakpoints inside unit te
"request": "launch",
"program": "${file}",
"purpose": ["debug-test"],
- "python": "/opt/miniconda3/envs/rmg_env/bin/python-jl",
+ "python": "/opt/miniconda3/envs/rmg_env/bin/python",
"console": "integratedTerminal",
"justMyCode": false,
"env": {"PYTEST_ADDOPTS": "--no-cov",} // without disabling coverage VS Code doesn't stop at breakpoints while debugging because pytest-cov is using the same technique to access the source code being run
@@ -255,7 +251,7 @@ To configure the rest of the settings, find the ``settings.json`` file in your `
You can use the following settings to configure the pytest framework::
"python.testing.pytestEnabled": true,
- "python.testing.pytestPath": "python-jl -m pytest",
+ "python.testing.pytestPath": "python -m pytest",
"python.testing.pytestArgs": [
"-p", "julia.pytestplugin",
"--julia-compiled-modules=no",
diff --git a/documentation/source/users/rmg/installation/dependencies.rst b/documentation/source/users/rmg/installation/dependencies.rst
index ba8a3ad66c..f91b280761 100644
--- a/documentation/source/users/rmg/installation/dependencies.rst
+++ b/documentation/source/users/rmg/installation/dependencies.rst
@@ -24,7 +24,6 @@ Briefly, RMG depends on the following packages, almost all of which can be found
* **graphviz:** generating flux diagrams
* **jinja2:** Python templating language for html rendering
* **jupyter:** (optional) for using IPython notebooks
-* **lpsolve:** mixed integer linear programming solver, used for resonance structure generation. Must also install Python extension.
* **markupsafe:** implements XML/HTML/XHTML markup safe strings for Python
* **matplotlib:** library for making plots
* **mock:** for unit-testing
diff --git a/documentation/source/users/rmg/running.rst b/documentation/source/users/rmg/running.rst
index 0bb21b80a6..aa28563a89 100755
--- a/documentation/source/users/rmg/running.rst
+++ b/documentation/source/users/rmg/running.rst
@@ -10,7 +10,7 @@ Running a basic RMG job is straightforward, as shown in the example below. Howev
Basic run::
- python-jl rmg.py input.py
+ python rmg.py input.py
.. _inputflags:
@@ -19,7 +19,7 @@ Input flags
The options for input flags can be found in ``/RMG-Py/rmgpy/util.py``. Running ::
- python-jl rmg.py -h
+ python rmg.py -h
at the command line will print the documentation from ``util.py``, which is reproduced below for convenience::
@@ -63,23 +63,23 @@ Some representative example usages are shown below.
Run by restarting from a seed mechanism::
- python-jl rmg.py -r path/to/seed/ input.py
+ python rmg.py -r path/to/seed/ input.py
Run with CPU time profiling::
- python-jl rmg.py -p input.py
+ python rmg.py -p input.py
Run with multiprocessing for reaction generation and QMTP::
- python-jl rmg.py -n input.py
+ python rmg.py -n input.py
Run with setting a limit on the maximum execution time::
- python-jl rmg.py -t input.py
+ python rmg.py -t input.py
Run with setting a limit on the maximum number of iterations::
- python-jl rmg.py -i input.py
+ python rmg.py -i input.py
Details on the multiprocessing implementation
diff --git a/environment.yml b/environment.yml
index da15b2dfe3..a5e272232d 100644
--- a/environment.yml
+++ b/environment.yml
@@ -18,92 +18,81 @@
# moved diffeqpy to pip and (temporarily) removed chemprop
# - April 17, 2024 Limit versions of cclib at advice of maintainers.
# - August 4, 2024 Restricted pyrms to <2
+# - May 14, 2024 Removed diffeqpy by switching to call SciMLBase and Sundials using JuliaCall
+# - March 15, 2024 - started migration to Python 3.9
+# - Mar 11, 2024 Removed Julia dependencies, now considered optional
+#
name: rmg_env
channels:
- - rmg
- conda-forge
- cantera
+ - rmg
dependencies:
# System-level dependencies - we could install these at the OS level
# but by installing them in the conda environment we get better control
- - cairo
- - cairocffi
- - ffmpeg
- - xlrd
- - xlwt
- - h5py
- - graphviz
- - markupsafe
- - psutil
+ - conda-forge::cairo
+ - conda-forge::cairocffi
+ - conda-forge::ffmpeg >= 7
+ - conda-forge::xlrd
+ - conda-forge::xlwt
+ - conda-forge::h5py
+ - conda-forge::graphviz >=12
+ - conda-forge::markupsafe
+ - conda-forge::psutil
# conda-forge not default, since default has a version information bug
# (see https://github.com/ReactionMechanismGenerator/RMG-Py/pull/2421)
- conda-forge::ncurses
- conda-forge::suitesparse
# external software tools for chemistry
- - coolprop
- - cantera::cantera=2.6
+ - conda-forge::coolprop
+ - cantera::cantera =2.6
- conda-forge::mopac
# see https://github.com/ReactionMechanismGenerator/RMG-Py/pull/2639#issuecomment-2050292972
- conda-forge::cclib >=1.6.3,<1.9
- conda-forge::openbabel >= 3
- conda-forge::rdkit >=2022.09.1
-# general-purpose external software tools
- - conda-forge::julia=1.9.1
- - conda-forge::pyjulia >=0.6
-
# Python tools
- - python >=3.7
- - coverage
- - cython >=0.25.2
- - scikit-learn
- - scipy <1.11
- - numpy >=1.10.0
- - pydot
- - jinja2
- - jupyter
- - pymongo
- - pyparsing
- - pyyaml
- - networkx
- - pytest
- - pytest-cov
- # we use a the pytest-check plugin, which is on Conda and PyPI, but the
- # version compatible with Python 3.7 is only on PyPI
- # switch to the conda version after upgrading to 3.11
- # - conda-forge::pytest-check
- - pip
- - pip:
- - pytest-check
- - matplotlib >=1.5
- - mpmath
- - pandas
+ - conda-forge::python >=3.9 # leave as GEQ so that GitHub actions can add EQ w/o breaking (contradictory deps)
+ - conda-forge::coverage
+ - conda-forge::cython >=0.25.2
+ - conda-forge::scikit-learn
+ - conda-forge::scipy >=1.9
+ - conda-forge::numpy >=1.10.0
+ - conda-forge::pydot
+ - conda-forge::jinja2
+ - conda-forge::jupyter
+ - conda-forge::pymongo
+ - conda-forge::pyparsing
+ - conda-forge::pyyaml
+ - conda-forge::networkx
+ - conda-forge::pytest
+ - conda-forge::pytest-cov
+ - conda-forge::pytest-check
+ - conda-forge::pyutilib
+ - conda-forge::matplotlib >=1.5
+ - conda-forge::mpmath
+ - conda-forge::pandas
- conda-forge::gprof2dot
- conda-forge::numdifftools
- - conda-forge::quantities
- - conda-forge::muq
- - conda-forge::lpsolve55
+ # bug in quantities, see:
+ # https://github.com/ReactionMechanismGenerator/RMG-Py/pull/2694#issuecomment-2489286263
+ - conda-forge::quantities !=0.16.0,!=0.16.1
- conda-forge::ringdecomposerlib-python
# packages we maintain
- rmg::pydas >=1.0.3
- rmg::pydqed >=1.0.3
- - rmg::pyrms <2
- rmg::symmetry
-# packages we would like to stop maintaining (and why)
- - rmg::diffeqpy
- # we should use the official verison https://github.com/SciML/diffeqpy),
- # rather than ours (which is only made so that we can get it from conda)
- # It is only on pip, so we will need to do something like:
- # https://stackoverflow.com/a/35245610
- # Note that _some other_ dep. in this list requires diffeqpy in its recipe
- # which will cause it to be downloaded from the rmg conda channel
-
# configure packages to use OpenBLAS instead of Intel MKL
- blas=*=openblas
+# optional dependencies for using ReactionMechanismSimulator
+# remove the leading '#' to install the required dependencies
+ # - conda-forge::pyjuliacall
+
# additional packages that are required, but not specified here (and why)
# pydqed, pydas, mopac, and likely others require a fortran compiler (specifically gfortran)
# in the environment. Normally we would add this to the environment file with
diff --git a/examples/arkane/networks/CH2NH2/input.py b/examples/arkane/networks/CH2NH2_mse/input.py
similarity index 100%
rename from examples/arkane/networks/CH2NH2/input.py
rename to examples/arkane/networks/CH2NH2_mse/input.py
diff --git a/examples/arkane/networks/CH2NH2/yaml_files/CH2NH.yml b/examples/arkane/networks/CH2NH2_mse/yaml_files/CH2NH.yml
similarity index 100%
rename from examples/arkane/networks/CH2NH2/yaml_files/CH2NH.yml
rename to examples/arkane/networks/CH2NH2_mse/yaml_files/CH2NH.yml
diff --git a/examples/arkane/networks/CH2NH2/yaml_files/CH2NH2.yml b/examples/arkane/networks/CH2NH2_mse/yaml_files/CH2NH2.yml
similarity index 100%
rename from examples/arkane/networks/CH2NH2/yaml_files/CH2NH2.yml
rename to examples/arkane/networks/CH2NH2_mse/yaml_files/CH2NH2.yml
diff --git a/examples/arkane/networks/CH2NH2/yaml_files/CH2NH_to_CH2NH2.yml b/examples/arkane/networks/CH2NH2_mse/yaml_files/CH2NH_to_CH2NH2.yml
similarity index 100%
rename from examples/arkane/networks/CH2NH2/yaml_files/CH2NH_to_CH2NH2.yml
rename to examples/arkane/networks/CH2NH2_mse/yaml_files/CH2NH_to_CH2NH2.yml
diff --git a/examples/arkane/networks/CH2NH2/yaml_files/CH2NH_to_CH3NH.yml b/examples/arkane/networks/CH2NH2_mse/yaml_files/CH2NH_to_CH3NH.yml
similarity index 100%
rename from examples/arkane/networks/CH2NH2/yaml_files/CH2NH_to_CH3NH.yml
rename to examples/arkane/networks/CH2NH2_mse/yaml_files/CH2NH_to_CH3NH.yml
diff --git a/examples/arkane/networks/CH2NH2/yaml_files/CH3NH.yml b/examples/arkane/networks/CH2NH2_mse/yaml_files/CH3NH.yml
similarity index 100%
rename from examples/arkane/networks/CH2NH2/yaml_files/CH3NH.yml
rename to examples/arkane/networks/CH2NH2_mse/yaml_files/CH3NH.yml
diff --git a/examples/arkane/networks/CH2NH2/yaml_files/CH3NH_to_CH2NH2.yml b/examples/arkane/networks/CH2NH2_mse/yaml_files/CH3NH_to_CH2NH2.yml
similarity index 100%
rename from examples/arkane/networks/CH2NH2/yaml_files/CH3NH_to_CH2NH2.yml
rename to examples/arkane/networks/CH2NH2_mse/yaml_files/CH3NH_to_CH2NH2.yml
diff --git a/examples/arkane/networks/CH2NH2/yaml_files/H.yml b/examples/arkane/networks/CH2NH2_mse/yaml_files/H.yml
similarity index 100%
rename from examples/arkane/networks/CH2NH2/yaml_files/H.yml
rename to examples/arkane/networks/CH2NH2_mse/yaml_files/H.yml
diff --git a/examples/arkane/networks/acetyl+O2/input.py b/examples/arkane/networks/acetyl+O2_cse/input.py
similarity index 100%
rename from examples/arkane/networks/acetyl+O2/input.py
rename to examples/arkane/networks/acetyl+O2_cse/input.py
diff --git a/examples/arkane/networks/acetyl+O2_rs/input.py b/examples/arkane/networks/acetyl+O2_rs/input.py
new file mode 100644
index 0000000000..84a57961fa
--- /dev/null
+++ b/examples/arkane/networks/acetyl+O2_rs/input.py
@@ -0,0 +1,284 @@
+################################################################################
+#
+# Arkane input file for acetyl + O2 pressure-dependent reaction network
+#
+################################################################################
+
+title = 'acetyl + oxygen'
+
+description = \
+"""
+The chemically-activated reaction of acetyl with oxygen. This system is of
+interest in atmospheric chemistry as a step in the conversion of acetaldehyde
+to the secondary pollutant peroxyacetylnitrate (PAN); it is also potentially
+important in the ignition chemistry of ethanol.
+"""
+
+species(
+ label = 'acetylperoxy',
+ structure = SMILES('CC(=O)O[O]'),
+ E0 = (-34.6,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(75.04,"g/mol")),
+ NonlinearRotor(inertia=([54.2977,104.836,156.05],"amu*angstrom^2"), symmetry=1),
+ HarmonicOscillator(frequencies=([319.695,500.474,536.674,543.894,727.156,973.365,1037.77,1119.72,1181.55,1391.11,1449.53,1454.72,1870.51,3037.12,3096.93,3136.39],"cm^-1")),
+ HinderedRotor(inertia=(7.38359,"amu*angstrom^2"), symmetry=1, fourier=([[-1.95191,-11.8215,0.740041,-0.049118,-0.464522],[0.000227764,0.00410782,-0.000805364,-0.000548218,-0.000266277]],"kJ/mol")),
+ HinderedRotor(inertia=(2.94723,"amu*angstrom^2"), symmetry=3, fourier=([[0.130647,0.0401507,-2.54582,-0.0436065,-0.120982],[-0.000701659,-0.000989654,0.00783349,-0.00140978,-0.00145843]],"kJ/mol")),
+ ],
+ spinMultiplicity = 2,
+ opticalIsomers = 1,
+ molecularWeight = (75.04,"g/mol"),
+ collisionModel = TransportData(sigma=(5.09,'angstrom'), epsilon=(473,'K')),
+ energyTransferModel = SingleExponentialDown(
+ alpha0 = (0.5718,'kcal/mol'),
+ T0 = (300,'K'),
+ n = 0.85,
+ ),
+)
+
+species(
+ label = 'hydroperoxylvinoxy',
+ structure = SMILES('[CH2]C(=O)OO'),
+ E0 = (-32.4,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(75.04,"g/mol")),
+ NonlinearRotor(inertia=([44.8034,110.225,155.029],"u*angstrom**2"), symmetry=1),
+ HarmonicOscillator(frequencies=([318.758,420.907,666.223,675.962,752.824,864.66,998.471,1019.54,1236.21,1437.91,1485.74,1687.9,3145.44,3262.88,3434.34],"cm^-1")),
+ HinderedRotor(inertia=(1.68464,"u*angstrom**2"), symmetry=2, fourier=([[0.359649,-16.1155,-0.593311,1.72918,0.256314],[-7.42981e-06,-0.000238057,3.29276e-05,-6.62608e-05,8.8443e-05]],"kJ/mol")),
+ HinderedRotor(inertia=(8.50433,"u*angstrom**2"), symmetry=1, fourier=([[-7.53504,-23.4471,-3.3974,-0.940593,-0.313674],[-4.58248,-2.47177,0.550012,1.03771,0.844226]],"kJ/mol")),
+ HinderedRotor(inertia=(0.803309,"u*angstrom**2"), symmetry=1, fourier=([[-8.65946,-3.97888,-1.13469,-0.402197,-0.145101],[4.41884e-05,4.83249e-05,1.30275e-05,-1.31353e-05,-6.66878e-06]],"kJ/mol")),
+ ],
+ spinMultiplicity = 2,
+ opticalIsomers = 1,
+ molecularWeight = (75.04,"g/mol"),
+ collisionModel = TransportData(sigma=(5.09,'angstrom'), epsilon=(473,'K')),
+ energyTransferModel = SingleExponentialDown(
+ alpha0 = (0.5718,'kcal/mol'),
+ T0 = (300,'K'),
+ n = 0.85,
+ ),
+)
+
+species(
+ label = 'acetyl',
+ structure = SMILES('C[C]=O'),
+ E0 = (0.0,'kcal/mol'), #(-20.5205,"kJ/mol")
+ modes = [
+ IdealGasTranslation(mass=(43.05,"g/mol")),
+ NonlinearRotor(inertia=([5.94518,50.8166,53.6436],"u*angstrom**2"), symmetry=1),
+ HarmonicOscillator(frequencies=([464.313,845.126,1010.54,1038.43,1343.54,1434.69,1442.25,1906.18,2985.46,3076.57,3079.46],"cm^-1")),
+ HinderedRotor(inertia=(1.61752,"u*angstrom**2"), symmetry=3, barrier=(2.00242,"kJ/mol")),
+ ],
+ spinMultiplicity = 2,
+ opticalIsomers = 1,
+)
+
+species(
+ label = 'oxygen',
+ structure = SMILES('[O][O]'),
+ E0 = (0.0,'kcal/mol'), #(-5.74557,"kJ/mol")
+ modes = [
+ IdealGasTranslation(mass=(32.00,"g/mol")),
+ LinearRotor(inertia=(11.6056,"u*angstrom**2"), symmetry=2),
+ HarmonicOscillator(frequencies=([1621.54],"cm^-1")),
+ ],
+ spinMultiplicity = 3,
+ opticalIsomers = 1,
+)
+
+species(
+ label = 'ketene',
+ structure = SMILES('C=C=O'),
+ E0 = (-6.6,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(42.0106,"g/mol")),
+ NonlinearRotor(inertia=([1.76922,48.8411,50.6103],"u*angstrom**2"), symmetry=2),
+ HarmonicOscillator(frequencies=([441.622,548.317,592.155,981.379,1159.66,1399.86,2192.1,3150.02,3240.58],"cm^-1")),
+ ],
+ spinMultiplicity = 1,
+ opticalIsomers = 1,
+)
+
+species(
+ label = 'lactone',
+ structure = SMILES('C1OC1(=O)'),
+ E0 = (-30.8,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(58.0055,"g/mol")),
+ NonlinearRotor(inertia=([20.1707,62.4381,79.1616],"u*angstrom**2"), symmetry=1),
+ HarmonicOscillator(frequencies=([484.387,527.771,705.332,933.372,985.563,1050.26,1107.93,1176.43,1466.54,1985.09,3086.23,3186.46],"cm^-1")),
+ ],
+ spinMultiplicity = 1,
+ opticalIsomers = 1,
+)
+
+
+
+species(
+ label = 'hydroxyl',
+ structure = SMILES('[OH]'),
+ E0 = (0.0,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(17.0027,"g/mol")),
+ LinearRotor(inertia=(0.899988,"u*angstrom**2"), symmetry=1),
+ HarmonicOscillator(frequencies=([3676.39],"cm^-1")),
+ ],
+ spinMultiplicity = 2,
+ opticalIsomers = 1,
+)
+
+species(
+ label = 'hydroperoxyl',
+ structure = SMILES('O[O]'),
+ E0 = (0.0,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(32.9977,"g/mol")),
+ NonlinearRotor(inertia=([0.811564,14.9434,15.755],"u*angstrom**2"), symmetry=1),
+ HarmonicOscillator(frequencies=([1156.77,1424.28,3571.81],"cm^-1")),
+ ],
+ spinMultiplicity = 2,
+ opticalIsomers = 1,
+)
+
+species(
+ label = 'nitrogen',
+ structure = SMILES('N#N'),
+ molecularWeight = (28.04,"g/mol"),
+ collisionModel = TransportData(sigma=(3.70,'angstrom'), epsilon=(94.9,'K')),
+ reactive = False
+)
+
+################################################################################
+
+transitionState(
+ label = 'entrance1',
+ E0 = (0.0,'kcal/mol'),
+)
+
+transitionState(
+ label = 'isom1',
+ E0 = (-5.8,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(75.04,"g/mol")),
+ NonlinearRotor(inertia=([49.3418,103.697,149.682],"u*angstrom**2"), symmetry=1, quantum=False),
+ HarmonicOscillator(frequencies=([148.551,306.791,484.573,536.709,599.366,675.538,832.594,918.413,1022.28,1031.45,1101.01,1130.05,1401.51,1701.26,1844.17,3078.6,3163.07],"cm^-1"), quantum=True),
+ ],
+ spinMultiplicity = 2,
+ opticalIsomers = 1,
+ frequency = (-1679.04,'cm^-1'),
+)
+
+transitionState(
+ label = 'exit1',
+ E0 = (0.6,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(75.04,"g/mol")),
+ NonlinearRotor(inertia=([55.4256, 136.1886, 188.2442],"u*angstrom**2"), symmetry=1),
+ HarmonicOscillator(frequencies=([59.9256,204.218,352.811,466.297,479.997,542.345,653.897,886.657,1017.91,1079.17,1250.02,1309.14,1370.36,1678.52,2162.41,3061.53,3135.55],"cm^-1")),
+ ],
+ spinMultiplicity = 2,
+ opticalIsomers = 1,
+ frequency=(-1053.25,'cm^-1'),
+)
+
+transitionState(
+ label = 'exit2',
+ E0 = (-4.6,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(75.0082,"g/mol"), quantum=False),
+ NonlinearRotor(inertia=([51.7432,119.373,171.117],"u*angstrom**2"), symmetry=1, quantum=False),
+ HarmonicOscillator(frequencies=([250.311,383.83,544.382,578.988,595.324,705.422,964.712,1103.25,1146.91,1415.98,1483.33,1983.79,3128,3143.84,3255.29],"cm^-1")),
+ HinderedRotor(inertia=(9.35921,"u*angstrom**2"), symmetry=1, fourier=([[-11.2387,-12.5928,-3.87844,1.13314,-0.358812],[-1.59863,-8.03329,-5.05235,3.13723,2.45989]],"kJ/mol")),
+ HinderedRotor(inertia=(0.754698,"u*angstrom**2"), symmetry=1, barrier=(47.7645,"kJ/mol")),
+ ],
+ spinMultiplicity = 2,
+ opticalIsomers = 1,
+ frequency=(-404.271,'cm^-1'),
+)
+
+transitionState(
+ label = 'exit3',
+ E0 = (-7.2,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(75.04,"g/mol")),
+ NonlinearRotor(inertia=([53.2821, 120.4050, 170.1570],"u*angstrom**2"), symmetry=1),
+ HarmonicOscillator(frequencies=([152.155,181.909,311.746,348.646,608.487,624.378,805.347,948.875,995.256,996.982,1169.1,1412.6,1834.83,3124.43,3245.2,3634.45],"cm^-1")),
+ HinderedRotor(inertia=(0.813269,"u*angstrom**2"), symmetry=1, fourier=([[-1.15338,-2.18664,-0.669531,-0.11502,-0.0512599],[0.00245222,0.0107485,0.00334564,-0.00165288,-0.0028674]],"kJ/mol")),
+ ],
+ spinMultiplicity = 2,
+ opticalIsomers = 1,
+ frequency=(-618.234,'cm^-1'),
+)
+
+################################################################################
+
+reaction(
+ label = 'entrance1',
+ reactants = ['acetyl', 'oxygen'],
+ products = ['acetylperoxy'],
+ transitionState = 'entrance1',
+ kinetics = Arrhenius(A=(2.65e6,'m^3/(mol*s)'), n=0.0, Ea=(0.0,'kcal/mol'), T0=(1,"K")),
+)
+
+reaction(
+ label = 'isom1',
+ reactants = ['acetylperoxy'],
+ products = ['hydroperoxylvinoxy'],
+ transitionState = 'isom1',
+ tunneling = 'Eckart',
+)
+
+reaction(
+ label = 'exit1',
+ reactants = ['acetylperoxy'],
+ products = ['ketene', 'hydroperoxyl'],
+ transitionState = 'exit1',
+ tunneling = 'Eckart',
+)
+
+reaction(
+ label = 'exit2',
+ reactants = ['hydroperoxylvinoxy'],
+ products = ['ketene', 'hydroperoxyl'],
+ transitionState = 'exit2',
+ tunneling = 'Eckart',
+)
+
+reaction(
+ label = 'exit3',
+ reactants = ['hydroperoxylvinoxy'],
+ products = ['lactone', 'hydroxyl'],
+ transitionState = 'exit3',
+ tunneling = 'Eckart',
+)
+
+#################################################################################
+
+network(
+ label = 'acetyl + O2',
+ isomers = [
+ 'acetylperoxy',
+ 'hydroperoxylvinoxy',
+ ],
+ reactants = [
+ ('acetyl', 'oxygen'),
+ ('ketene', 'hydroperoxyl'),
+ ('lactone', 'hydroxyl')
+ ],
+ bathGas = {
+ 'nitrogen': 1.0,
+ }
+)
+
+#################################################################################
+
+pressureDependence(
+ 'acetyl + O2',
+ Tmin=(300.0,'K'), Tmax=(2000.0,'K'), Tcount=8,
+ Pmin=(0.01,'bar'), Pmax=(100.0,'bar'), Pcount=5,
+ maximumGrainSize = (1.0,'kcal/mol'),
+ minimumGrainCount = 250,
+ method = 'reservoir state',
+ interpolationModel = ('chebyshev', 6, 4),
+ activeJRotor = True,
+)
diff --git a/examples/arkane/networks/acetyl+O2_sls/input.py b/examples/arkane/networks/acetyl+O2_sls/input.py
new file mode 100644
index 0000000000..4595e288fc
--- /dev/null
+++ b/examples/arkane/networks/acetyl+O2_sls/input.py
@@ -0,0 +1,284 @@
+################################################################################
+#
+# Arkane input file for acetyl + O2 pressure-dependent reaction network
+#
+################################################################################
+
+title = 'acetyl + oxygen'
+
+description = \
+"""
+The chemically-activated reaction of acetyl with oxygen. This system is of
+interest in atmospheric chemistry as a step in the conversion of acetaldehyde
+to the secondary pollutant peroxyacetylnitrate (PAN); it is also potentially
+important in the ignition chemistry of ethanol.
+"""
+
+species(
+ label = 'acetylperoxy',
+ structure = SMILES('CC(=O)O[O]'),
+ E0 = (-34.6,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(75.04,"g/mol")),
+ NonlinearRotor(inertia=([54.2977,104.836,156.05],"amu*angstrom^2"), symmetry=1),
+ HarmonicOscillator(frequencies=([319.695,500.474,536.674,543.894,727.156,973.365,1037.77,1119.72,1181.55,1391.11,1449.53,1454.72,1870.51,3037.12,3096.93,3136.39],"cm^-1")),
+ HinderedRotor(inertia=(7.38359,"amu*angstrom^2"), symmetry=1, fourier=([[-1.95191,-11.8215,0.740041,-0.049118,-0.464522],[0.000227764,0.00410782,-0.000805364,-0.000548218,-0.000266277]],"kJ/mol")),
+ HinderedRotor(inertia=(2.94723,"amu*angstrom^2"), symmetry=3, fourier=([[0.130647,0.0401507,-2.54582,-0.0436065,-0.120982],[-0.000701659,-0.000989654,0.00783349,-0.00140978,-0.00145843]],"kJ/mol")),
+ ],
+ spinMultiplicity = 2,
+ opticalIsomers = 1,
+ molecularWeight = (75.04,"g/mol"),
+ collisionModel = TransportData(sigma=(5.09,'angstrom'), epsilon=(473,'K')),
+ energyTransferModel = SingleExponentialDown(
+ alpha0 = (0.5718,'kcal/mol'),
+ T0 = (300,'K'),
+ n = 0.85,
+ ),
+)
+
+species(
+ label = 'hydroperoxylvinoxy',
+ structure = SMILES('[CH2]C(=O)OO'),
+ E0 = (-32.4,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(75.04,"g/mol")),
+ NonlinearRotor(inertia=([44.8034,110.225,155.029],"u*angstrom**2"), symmetry=1),
+ HarmonicOscillator(frequencies=([318.758,420.907,666.223,675.962,752.824,864.66,998.471,1019.54,1236.21,1437.91,1485.74,1687.9,3145.44,3262.88,3434.34],"cm^-1")),
+ HinderedRotor(inertia=(1.68464,"u*angstrom**2"), symmetry=2, fourier=([[0.359649,-16.1155,-0.593311,1.72918,0.256314],[-7.42981e-06,-0.000238057,3.29276e-05,-6.62608e-05,8.8443e-05]],"kJ/mol")),
+ HinderedRotor(inertia=(8.50433,"u*angstrom**2"), symmetry=1, fourier=([[-7.53504,-23.4471,-3.3974,-0.940593,-0.313674],[-4.58248,-2.47177,0.550012,1.03771,0.844226]],"kJ/mol")),
+ HinderedRotor(inertia=(0.803309,"u*angstrom**2"), symmetry=1, fourier=([[-8.65946,-3.97888,-1.13469,-0.402197,-0.145101],[4.41884e-05,4.83249e-05,1.30275e-05,-1.31353e-05,-6.66878e-06]],"kJ/mol")),
+ ],
+ spinMultiplicity = 2,
+ opticalIsomers = 1,
+ molecularWeight = (75.04,"g/mol"),
+ collisionModel = TransportData(sigma=(5.09,'angstrom'), epsilon=(473,'K')),
+ energyTransferModel = SingleExponentialDown(
+ alpha0 = (0.5718,'kcal/mol'),
+ T0 = (300,'K'),
+ n = 0.85,
+ ),
+)
+
+species(
+ label = 'acetyl',
+ structure = SMILES('C[C]=O'),
+ E0 = (0.0,'kcal/mol'), #(-20.5205,"kJ/mol")
+ modes = [
+ IdealGasTranslation(mass=(43.05,"g/mol")),
+ NonlinearRotor(inertia=([5.94518,50.8166,53.6436],"u*angstrom**2"), symmetry=1),
+ HarmonicOscillator(frequencies=([464.313,845.126,1010.54,1038.43,1343.54,1434.69,1442.25,1906.18,2985.46,3076.57,3079.46],"cm^-1")),
+ HinderedRotor(inertia=(1.61752,"u*angstrom**2"), symmetry=3, barrier=(2.00242,"kJ/mol")),
+ ],
+ spinMultiplicity = 2,
+ opticalIsomers = 1,
+)
+
+species(
+ label = 'oxygen',
+ structure = SMILES('[O][O]'),
+ E0 = (0.0,'kcal/mol'), #(-5.74557,"kJ/mol")
+ modes = [
+ IdealGasTranslation(mass=(32.00,"g/mol")),
+ LinearRotor(inertia=(11.6056,"u*angstrom**2"), symmetry=2),
+ HarmonicOscillator(frequencies=([1621.54],"cm^-1")),
+ ],
+ spinMultiplicity = 3,
+ opticalIsomers = 1,
+)
+
+species(
+ label = 'ketene',
+ structure = SMILES('C=C=O'),
+ E0 = (-6.6,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(42.0106,"g/mol")),
+ NonlinearRotor(inertia=([1.76922,48.8411,50.6103],"u*angstrom**2"), symmetry=2),
+ HarmonicOscillator(frequencies=([441.622,548.317,592.155,981.379,1159.66,1399.86,2192.1,3150.02,3240.58],"cm^-1")),
+ ],
+ spinMultiplicity = 1,
+ opticalIsomers = 1,
+)
+
+species(
+ label = 'lactone',
+ structure = SMILES('C1OC1(=O)'),
+ E0 = (-30.8,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(58.0055,"g/mol")),
+ NonlinearRotor(inertia=([20.1707,62.4381,79.1616],"u*angstrom**2"), symmetry=1),
+ HarmonicOscillator(frequencies=([484.387,527.771,705.332,933.372,985.563,1050.26,1107.93,1176.43,1466.54,1985.09,3086.23,3186.46],"cm^-1")),
+ ],
+ spinMultiplicity = 1,
+ opticalIsomers = 1,
+)
+
+
+
+species(
+ label = 'hydroxyl',
+ structure = SMILES('[OH]'),
+ E0 = (0.0,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(17.0027,"g/mol")),
+ LinearRotor(inertia=(0.899988,"u*angstrom**2"), symmetry=1),
+ HarmonicOscillator(frequencies=([3676.39],"cm^-1")),
+ ],
+ spinMultiplicity = 2,
+ opticalIsomers = 1,
+)
+
+species(
+ label = 'hydroperoxyl',
+ structure = SMILES('O[O]'),
+ E0 = (0.0,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(32.9977,"g/mol")),
+ NonlinearRotor(inertia=([0.811564,14.9434,15.755],"u*angstrom**2"), symmetry=1),
+ HarmonicOscillator(frequencies=([1156.77,1424.28,3571.81],"cm^-1")),
+ ],
+ spinMultiplicity = 2,
+ opticalIsomers = 1,
+)
+
+species(
+ label = 'nitrogen',
+ structure = SMILES('N#N'),
+ molecularWeight = (28.04,"g/mol"),
+ collisionModel = TransportData(sigma=(3.70,'angstrom'), epsilon=(94.9,'K')),
+ reactive = False
+)
+
+################################################################################
+
+transitionState(
+ label = 'entrance1',
+ E0 = (0.0,'kcal/mol'),
+)
+
+transitionState(
+ label = 'isom1',
+ E0 = (-5.8,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(75.04,"g/mol")),
+ NonlinearRotor(inertia=([49.3418,103.697,149.682],"u*angstrom**2"), symmetry=1, quantum=False),
+ HarmonicOscillator(frequencies=([148.551,306.791,484.573,536.709,599.366,675.538,832.594,918.413,1022.28,1031.45,1101.01,1130.05,1401.51,1701.26,1844.17,3078.6,3163.07],"cm^-1"), quantum=True),
+ ],
+ spinMultiplicity = 2,
+ opticalIsomers = 1,
+ frequency = (-1679.04,'cm^-1'),
+)
+
+transitionState(
+ label = 'exit1',
+ E0 = (0.6,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(75.04,"g/mol")),
+ NonlinearRotor(inertia=([55.4256, 136.1886, 188.2442],"u*angstrom**2"), symmetry=1),
+ HarmonicOscillator(frequencies=([59.9256,204.218,352.811,466.297,479.997,542.345,653.897,886.657,1017.91,1079.17,1250.02,1309.14,1370.36,1678.52,2162.41,3061.53,3135.55],"cm^-1")),
+ ],
+ spinMultiplicity = 2,
+ opticalIsomers = 1,
+ frequency=(-1053.25,'cm^-1'),
+)
+
+transitionState(
+ label = 'exit2',
+ E0 = (-4.6,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(75.0082,"g/mol"), quantum=False),
+ NonlinearRotor(inertia=([51.7432,119.373,171.117],"u*angstrom**2"), symmetry=1, quantum=False),
+ HarmonicOscillator(frequencies=([250.311,383.83,544.382,578.988,595.324,705.422,964.712,1103.25,1146.91,1415.98,1483.33,1983.79,3128,3143.84,3255.29],"cm^-1")),
+ HinderedRotor(inertia=(9.35921,"u*angstrom**2"), symmetry=1, fourier=([[-11.2387,-12.5928,-3.87844,1.13314,-0.358812],[-1.59863,-8.03329,-5.05235,3.13723,2.45989]],"kJ/mol")),
+ HinderedRotor(inertia=(0.754698,"u*angstrom**2"), symmetry=1, barrier=(47.7645,"kJ/mol")),
+ ],
+ spinMultiplicity = 2,
+ opticalIsomers = 1,
+ frequency=(-404.271,'cm^-1'),
+)
+
+transitionState(
+ label = 'exit3',
+ E0 = (-7.2,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(75.04,"g/mol")),
+ NonlinearRotor(inertia=([53.2821, 120.4050, 170.1570],"u*angstrom**2"), symmetry=1),
+ HarmonicOscillator(frequencies=([152.155,181.909,311.746,348.646,608.487,624.378,805.347,948.875,995.256,996.982,1169.1,1412.6,1834.83,3124.43,3245.2,3634.45],"cm^-1")),
+ HinderedRotor(inertia=(0.813269,"u*angstrom**2"), symmetry=1, fourier=([[-1.15338,-2.18664,-0.669531,-0.11502,-0.0512599],[0.00245222,0.0107485,0.00334564,-0.00165288,-0.0028674]],"kJ/mol")),
+ ],
+ spinMultiplicity = 2,
+ opticalIsomers = 1,
+ frequency=(-618.234,'cm^-1'),
+)
+
+################################################################################
+
+reaction(
+ label = 'entrance1',
+ reactants = ['acetyl', 'oxygen'],
+ products = ['acetylperoxy'],
+ transitionState = 'entrance1',
+ kinetics = Arrhenius(A=(2.65e6,'m^3/(mol*s)'), n=0.0, Ea=(0.0,'kcal/mol'), T0=(1,"K")),
+)
+
+reaction(
+ label = 'isom1',
+ reactants = ['acetylperoxy'],
+ products = ['hydroperoxylvinoxy'],
+ transitionState = 'isom1',
+ tunneling = 'Eckart',
+)
+
+reaction(
+ label = 'exit1',
+ reactants = ['acetylperoxy'],
+ products = ['ketene', 'hydroperoxyl'],
+ transitionState = 'exit1',
+ tunneling = 'Eckart',
+)
+
+reaction(
+ label = 'exit2',
+ reactants = ['hydroperoxylvinoxy'],
+ products = ['ketene', 'hydroperoxyl'],
+ transitionState = 'exit2',
+ tunneling = 'Eckart',
+)
+
+reaction(
+ label = 'exit3',
+ reactants = ['hydroperoxylvinoxy'],
+ products = ['lactone', 'hydroxyl'],
+ transitionState = 'exit3',
+ tunneling = 'Eckart',
+)
+
+#################################################################################
+
+network(
+ label = 'acetyl + O2',
+ isomers = [
+ 'acetylperoxy',
+ 'hydroperoxylvinoxy',
+ ],
+ reactants = [
+ ('acetyl', 'oxygen'),
+ ('ketene', 'hydroperoxyl'),
+ ('lactone', 'hydroxyl')
+ ],
+ bathGas = {
+ 'nitrogen': 1.0,
+ }
+)
+
+#################################################################################
+
+pressureDependence(
+ 'acetyl + O2',
+ Tmin=(300.0,'K'), Tmax=(2000.0,'K'), Tcount=8,
+ Pmin=(0.01,'bar'), Pmax=(100.0,'bar'), Pcount=5,
+ maximumGrainSize = (1.0,'kcal/mol'),
+ minimumGrainCount = 250,
+ method = 'simulation least squares',
+ interpolationModel = ('chebyshev', 6, 4),
+ activeJRotor = True,
+)
diff --git a/examples/arkane/networks/acetyl+O2_slse/input.py b/examples/arkane/networks/acetyl+O2_slse/input.py
new file mode 100644
index 0000000000..46114439b4
--- /dev/null
+++ b/examples/arkane/networks/acetyl+O2_slse/input.py
@@ -0,0 +1,284 @@
+################################################################################
+#
+# Arkane input file for acetyl + O2 pressure-dependent reaction network
+#
+################################################################################
+
+title = 'acetyl + oxygen'
+
+description = \
+"""
+The chemically-activated reaction of acetyl with oxygen. This system is of
+interest in atmospheric chemistry as a step in the conversion of acetaldehyde
+to the secondary pollutant peroxyacetylnitrate (PAN); it is also potentially
+important in the ignition chemistry of ethanol.
+"""
+
+species(
+ label = 'acetylperoxy',
+ structure = SMILES('CC(=O)O[O]'),
+ E0 = (-34.6,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(75.04,"g/mol")),
+ NonlinearRotor(inertia=([54.2977,104.836,156.05],"amu*angstrom^2"), symmetry=1),
+ HarmonicOscillator(frequencies=([319.695,500.474,536.674,543.894,727.156,973.365,1037.77,1119.72,1181.55,1391.11,1449.53,1454.72,1870.51,3037.12,3096.93,3136.39],"cm^-1")),
+ HinderedRotor(inertia=(7.38359,"amu*angstrom^2"), symmetry=1, fourier=([[-1.95191,-11.8215,0.740041,-0.049118,-0.464522],[0.000227764,0.00410782,-0.000805364,-0.000548218,-0.000266277]],"kJ/mol")),
+ HinderedRotor(inertia=(2.94723,"amu*angstrom^2"), symmetry=3, fourier=([[0.130647,0.0401507,-2.54582,-0.0436065,-0.120982],[-0.000701659,-0.000989654,0.00783349,-0.00140978,-0.00145843]],"kJ/mol")),
+ ],
+ spinMultiplicity = 2,
+ opticalIsomers = 1,
+ molecularWeight = (75.04,"g/mol"),
+ collisionModel = TransportData(sigma=(5.09,'angstrom'), epsilon=(473,'K')),
+ energyTransferModel = SingleExponentialDown(
+ alpha0 = (0.5718,'kcal/mol'),
+ T0 = (300,'K'),
+ n = 0.85,
+ ),
+)
+
+species(
+ label = 'hydroperoxylvinoxy',
+ structure = SMILES('[CH2]C(=O)OO'),
+ E0 = (-32.4,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(75.04,"g/mol")),
+ NonlinearRotor(inertia=([44.8034,110.225,155.029],"u*angstrom**2"), symmetry=1),
+ HarmonicOscillator(frequencies=([318.758,420.907,666.223,675.962,752.824,864.66,998.471,1019.54,1236.21,1437.91,1485.74,1687.9,3145.44,3262.88,3434.34],"cm^-1")),
+ HinderedRotor(inertia=(1.68464,"u*angstrom**2"), symmetry=2, fourier=([[0.359649,-16.1155,-0.593311,1.72918,0.256314],[-7.42981e-06,-0.000238057,3.29276e-05,-6.62608e-05,8.8443e-05]],"kJ/mol")),
+ HinderedRotor(inertia=(8.50433,"u*angstrom**2"), symmetry=1, fourier=([[-7.53504,-23.4471,-3.3974,-0.940593,-0.313674],[-4.58248,-2.47177,0.550012,1.03771,0.844226]],"kJ/mol")),
+ HinderedRotor(inertia=(0.803309,"u*angstrom**2"), symmetry=1, fourier=([[-8.65946,-3.97888,-1.13469,-0.402197,-0.145101],[4.41884e-05,4.83249e-05,1.30275e-05,-1.31353e-05,-6.66878e-06]],"kJ/mol")),
+ ],
+ spinMultiplicity = 2,
+ opticalIsomers = 1,
+ molecularWeight = (75.04,"g/mol"),
+ collisionModel = TransportData(sigma=(5.09,'angstrom'), epsilon=(473,'K')),
+ energyTransferModel = SingleExponentialDown(
+ alpha0 = (0.5718,'kcal/mol'),
+ T0 = (300,'K'),
+ n = 0.85,
+ ),
+)
+
+species(
+ label = 'acetyl',
+ structure = SMILES('C[C]=O'),
+ E0 = (0.0,'kcal/mol'), #(-20.5205,"kJ/mol")
+ modes = [
+ IdealGasTranslation(mass=(43.05,"g/mol")),
+ NonlinearRotor(inertia=([5.94518,50.8166,53.6436],"u*angstrom**2"), symmetry=1),
+ HarmonicOscillator(frequencies=([464.313,845.126,1010.54,1038.43,1343.54,1434.69,1442.25,1906.18,2985.46,3076.57,3079.46],"cm^-1")),
+ HinderedRotor(inertia=(1.61752,"u*angstrom**2"), symmetry=3, barrier=(2.00242,"kJ/mol")),
+ ],
+ spinMultiplicity = 2,
+ opticalIsomers = 1,
+)
+
+species(
+ label = 'oxygen',
+ structure = SMILES('[O][O]'),
+ E0 = (0.0,'kcal/mol'), #(-5.74557,"kJ/mol")
+ modes = [
+ IdealGasTranslation(mass=(32.00,"g/mol")),
+ LinearRotor(inertia=(11.6056,"u*angstrom**2"), symmetry=2),
+ HarmonicOscillator(frequencies=([1621.54],"cm^-1")),
+ ],
+ spinMultiplicity = 3,
+ opticalIsomers = 1,
+)
+
+species(
+ label = 'ketene',
+ structure = SMILES('C=C=O'),
+ E0 = (-6.6,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(42.0106,"g/mol")),
+ NonlinearRotor(inertia=([1.76922,48.8411,50.6103],"u*angstrom**2"), symmetry=2),
+ HarmonicOscillator(frequencies=([441.622,548.317,592.155,981.379,1159.66,1399.86,2192.1,3150.02,3240.58],"cm^-1")),
+ ],
+ spinMultiplicity = 1,
+ opticalIsomers = 1,
+)
+
+species(
+ label = 'lactone',
+ structure = SMILES('C1OC1(=O)'),
+ E0 = (-30.8,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(58.0055,"g/mol")),
+ NonlinearRotor(inertia=([20.1707,62.4381,79.1616],"u*angstrom**2"), symmetry=1),
+ HarmonicOscillator(frequencies=([484.387,527.771,705.332,933.372,985.563,1050.26,1107.93,1176.43,1466.54,1985.09,3086.23,3186.46],"cm^-1")),
+ ],
+ spinMultiplicity = 1,
+ opticalIsomers = 1,
+)
+
+
+
+species(
+ label = 'hydroxyl',
+ structure = SMILES('[OH]'),
+ E0 = (0.0,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(17.0027,"g/mol")),
+ LinearRotor(inertia=(0.899988,"u*angstrom**2"), symmetry=1),
+ HarmonicOscillator(frequencies=([3676.39],"cm^-1")),
+ ],
+ spinMultiplicity = 2,
+ opticalIsomers = 1,
+)
+
+species(
+ label = 'hydroperoxyl',
+ structure = SMILES('O[O]'),
+ E0 = (0.0,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(32.9977,"g/mol")),
+ NonlinearRotor(inertia=([0.811564,14.9434,15.755],"u*angstrom**2"), symmetry=1),
+ HarmonicOscillator(frequencies=([1156.77,1424.28,3571.81],"cm^-1")),
+ ],
+ spinMultiplicity = 2,
+ opticalIsomers = 1,
+)
+
+species(
+ label = 'nitrogen',
+ structure = SMILES('N#N'),
+ molecularWeight = (28.04,"g/mol"),
+ collisionModel = TransportData(sigma=(3.70,'angstrom'), epsilon=(94.9,'K')),
+ reactive = False
+)
+
+################################################################################
+
+transitionState(
+ label = 'entrance1',
+ E0 = (0.0,'kcal/mol'),
+)
+
+transitionState(
+ label = 'isom1',
+ E0 = (-5.8,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(75.04,"g/mol")),
+ NonlinearRotor(inertia=([49.3418,103.697,149.682],"u*angstrom**2"), symmetry=1, quantum=False),
+ HarmonicOscillator(frequencies=([148.551,306.791,484.573,536.709,599.366,675.538,832.594,918.413,1022.28,1031.45,1101.01,1130.05,1401.51,1701.26,1844.17,3078.6,3163.07],"cm^-1"), quantum=True),
+ ],
+ spinMultiplicity = 2,
+ opticalIsomers = 1,
+ frequency = (-1679.04,'cm^-1'),
+)
+
+transitionState(
+ label = 'exit1',
+ E0 = (0.6,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(75.04,"g/mol")),
+ NonlinearRotor(inertia=([55.4256, 136.1886, 188.2442],"u*angstrom**2"), symmetry=1),
+ HarmonicOscillator(frequencies=([59.9256,204.218,352.811,466.297,479.997,542.345,653.897,886.657,1017.91,1079.17,1250.02,1309.14,1370.36,1678.52,2162.41,3061.53,3135.55],"cm^-1")),
+ ],
+ spinMultiplicity = 2,
+ opticalIsomers = 1,
+ frequency=(-1053.25,'cm^-1'),
+)
+
+transitionState(
+ label = 'exit2',
+ E0 = (-4.6,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(75.0082,"g/mol"), quantum=False),
+ NonlinearRotor(inertia=([51.7432,119.373,171.117],"u*angstrom**2"), symmetry=1, quantum=False),
+ HarmonicOscillator(frequencies=([250.311,383.83,544.382,578.988,595.324,705.422,964.712,1103.25,1146.91,1415.98,1483.33,1983.79,3128,3143.84,3255.29],"cm^-1")),
+ HinderedRotor(inertia=(9.35921,"u*angstrom**2"), symmetry=1, fourier=([[-11.2387,-12.5928,-3.87844,1.13314,-0.358812],[-1.59863,-8.03329,-5.05235,3.13723,2.45989]],"kJ/mol")),
+ HinderedRotor(inertia=(0.754698,"u*angstrom**2"), symmetry=1, barrier=(47.7645,"kJ/mol")),
+ ],
+ spinMultiplicity = 2,
+ opticalIsomers = 1,
+ frequency=(-404.271,'cm^-1'),
+)
+
+transitionState(
+ label = 'exit3',
+ E0 = (-7.2,'kcal/mol'),
+ modes = [
+ IdealGasTranslation(mass=(75.04,"g/mol")),
+ NonlinearRotor(inertia=([53.2821, 120.4050, 170.1570],"u*angstrom**2"), symmetry=1),
+ HarmonicOscillator(frequencies=([152.155,181.909,311.746,348.646,608.487,624.378,805.347,948.875,995.256,996.982,1169.1,1412.6,1834.83,3124.43,3245.2,3634.45],"cm^-1")),
+ HinderedRotor(inertia=(0.813269,"u*angstrom**2"), symmetry=1, fourier=([[-1.15338,-2.18664,-0.669531,-0.11502,-0.0512599],[0.00245222,0.0107485,0.00334564,-0.00165288,-0.0028674]],"kJ/mol")),
+ ],
+ spinMultiplicity = 2,
+ opticalIsomers = 1,
+ frequency=(-618.234,'cm^-1'),
+)
+
+################################################################################
+
+reaction(
+ label = 'entrance1',
+ reactants = ['acetyl', 'oxygen'],
+ products = ['acetylperoxy'],
+ transitionState = 'entrance1',
+ kinetics = Arrhenius(A=(2.65e6,'m^3/(mol*s)'), n=0.0, Ea=(0.0,'kcal/mol'), T0=(1,"K")),
+)
+
+reaction(
+ label = 'isom1',
+ reactants = ['acetylperoxy'],
+ products = ['hydroperoxylvinoxy'],
+ transitionState = 'isom1',
+ tunneling = 'Eckart',
+)
+
+reaction(
+ label = 'exit1',
+ reactants = ['acetylperoxy'],
+ products = ['ketene', 'hydroperoxyl'],
+ transitionState = 'exit1',
+ tunneling = 'Eckart',
+)
+
+reaction(
+ label = 'exit2',
+ reactants = ['hydroperoxylvinoxy'],
+ products = ['ketene', 'hydroperoxyl'],
+ transitionState = 'exit2',
+ tunneling = 'Eckart',
+)
+
+reaction(
+ label = 'exit3',
+ reactants = ['hydroperoxylvinoxy'],
+ products = ['lactone', 'hydroxyl'],
+ transitionState = 'exit3',
+ tunneling = 'Eckart',
+)
+
+#################################################################################
+
+network(
+ label = 'acetyl + O2',
+ isomers = [
+ 'acetylperoxy',
+ 'hydroperoxylvinoxy',
+ ],
+ reactants = [
+ ('acetyl', 'oxygen'),
+ ('ketene', 'hydroperoxyl'),
+ ('lactone', 'hydroxyl')
+ ],
+ bathGas = {
+ 'nitrogen': 1.0,
+ }
+)
+
+#################################################################################
+
+pressureDependence(
+ 'acetyl + O2',
+ Tmin=(300.0,'K'), Tmax=(2000.0,'K'), Tcount=8,
+ Pmin=(0.01,'bar'), Pmax=(100.0,'bar'), Pcount=5,
+ maximumGrainSize = (1.0,'kcal/mol'),
+ minimumGrainCount = 250,
+ method = 'simulation least squares eigen',
+ interpolationModel = ('chebyshev', 6, 4),
+ activeJRotor = True,
+)
diff --git a/examples/arkane/networks/n-butanol/input.py b/examples/arkane/networks/n-butanol_msc/input.py
similarity index 100%
rename from examples/arkane/networks/n-butanol/input.py
rename to examples/arkane/networks/n-butanol_msc/input.py
diff --git a/examples/rmg/catalysis/ch4_o2/simulate_ch4o2cat_RMS.ipynb b/examples/rmg/catalysis/ch4_o2/simulate_ch4o2cat_RMS.ipynb
index 700eb3cac3..371bfc2cad 100644
--- a/examples/rmg/catalysis/ch4_o2/simulate_ch4o2cat_RMS.ipynb
+++ b/examples/rmg/catalysis/ch4_o2/simulate_ch4o2cat_RMS.ipynb
@@ -49,7 +49,7 @@
"metadata": {},
"outputs": [],
"source": [
- "# ! python-jl /rmg/RMG-Py/rmg.py /rmg/RMG-Py/examples/rmg/catalysis/ch4_o2/input.py"
+ "# ! python /rmg/RMG-Py/rmg.py /rmg/RMG-Py/examples/rmg/catalysis/ch4_o2/input.py"
]
},
{
diff --git a/examples/rmg/superminimal/sensitivity_analysis_superminimal_RMS.ipynb b/examples/rmg/superminimal/sensitivity_analysis_superminimal_RMS.ipynb
index ba6b3b6a98..44befc780f 100644
--- a/examples/rmg/superminimal/sensitivity_analysis_superminimal_RMS.ipynb
+++ b/examples/rmg/superminimal/sensitivity_analysis_superminimal_RMS.ipynb
@@ -40,7 +40,7 @@
"metadata": {},
"outputs": [],
"source": [
- "# ! python-jl /rmg/RMG-Py/rmg.py /rmg/RMG-Py/examples/rmg/superminimal/input.py"
+ "# ! python /rmg/RMG-Py/rmg.py /rmg/RMG-Py/examples/rmg/superminimal/input.py"
]
},
{
diff --git a/examples/rmg/superminimal/simulation_flux_rop_superminimal_RMS.ipynb b/examples/rmg/superminimal/simulation_flux_rop_superminimal_RMS.ipynb
index 7540b21b01..8ca7eb0029 100644
--- a/examples/rmg/superminimal/simulation_flux_rop_superminimal_RMS.ipynb
+++ b/examples/rmg/superminimal/simulation_flux_rop_superminimal_RMS.ipynb
@@ -40,7 +40,7 @@
"metadata": {},
"outputs": [],
"source": [
- "# ! python-jl /rmg/RMG-Py/rmg.py /rmg/RMG-Py/examples/rmg/superminimal/input.py"
+ "# ! python /rmg/RMG-Py/rmg.py /rmg/RMG-Py/examples/rmg/superminimal/input.py"
]
},
{
diff --git a/pytest.ini b/pytest.ini
index b17574cc17..f7643013e8 100644
--- a/pytest.ini
+++ b/pytest.ini
@@ -17,7 +17,7 @@ filterwarnings =
# -n auto # will use all available cores, if you first pip install pytest-xdist, but it doesn't work well with julia (so add --no-julia? or give up?)
addopts =
--keep-duplicates
- -p julia.pytestplugin
+ -s
-vv
--ignore test/regression
--cov=arkane --cov=rmgpy
diff --git a/rmg.py b/rmg.py
index 26c27c00a4..0b27098b81 100644
--- a/rmg.py
+++ b/rmg.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python-jl
+#!/usr/bin/env python
###############################################################################
# #
diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py
index 43e3836169..ab867f71d7 100644
--- a/rmgpy/data/kinetics/family.py
+++ b/rmgpy/data/kinetics/family.py
@@ -3535,8 +3535,8 @@ def make_bm_rules_from_template_rxn_map(self, template_rxn_map, nprocs=1, Tref=1
entries = list(self.groups.entries.values())
rxnlists = [(template_rxn_map[entry.label], entry.label)
if entry.label in template_rxn_map.keys() else [] for entry in entries]
- inputs = np.array([(self.forward_recipe.actions, rxns, Tref, fmax, label, [r.rank for r in rxns])
- for rxns, label in rxnlists])
+ inputs = [(self.forward_recipe.actions, rxns, Tref, fmax, label, [r.rank for r in rxns])
+ for rxns, label in rxnlists]
inds = np.arange(len(inputs))
np.random.shuffle(inds) # want to parallelize in random order
@@ -3545,9 +3545,9 @@ def make_bm_rules_from_template_rxn_map(self, template_rxn_map, nprocs=1, Tref=1
if nprocs > 1:
pool = mp.Pool(nprocs)
- kinetics_list = np.array(pool.map(_make_rule, inputs[inds]))
+ kinetics_list = np.array(pool.map(_make_rule, list(inputs[i] for i in inds)))
else:
- kinetics_list = np.array(list(map(_make_rule, inputs[inds])))
+ kinetics_list = np.array(list(map(_make_rule, list(inputs[i] for i in inds))))
kinetics_list = kinetics_list[revinds] # fix order
diff --git a/rmgpy/exceptions.py b/rmgpy/exceptions.py
index 7b4420cf66..f2e7dae460 100644
--- a/rmgpy/exceptions.py
+++ b/rmgpy/exceptions.py
@@ -109,15 +109,6 @@ class ForbiddenStructureException(Exception):
pass
-class ILPSolutionError(Exception):
- """
- An exception to be raised when solving an integer linear programming problem if a solution
- could not be found or the solution is not valid. Can pass a string to indicate the reason
- that the solution is invalid.
- """
- pass
-
-
class ImplicitBenzeneError(Exception):
"""
An exception class when encountering a group with too many implicit benzene
diff --git a/rmgpy/kinetics/kineticsdata.pyx b/rmgpy/kinetics/kineticsdata.pyx
index 7b4617cdd0..292a56a7ab 100644
--- a/rmgpy/kinetics/kineticsdata.pyx
+++ b/rmgpy/kinetics/kineticsdata.pyx
@@ -249,9 +249,8 @@ cdef class PDepKineticsData(PDepKineticsModel):
Plow = Pdata[j]
Phigh = Pdata[j + 1]
if Plow <= P and P <= Phigh:
- klow = kdata[i, j] * (kdata[i + 1, j] / kdata[i, j]) ** ((T - Tlow) / (Thigh - Tlow))
- khigh = kdata[i, j + 1] * (kdata[i + 1, j + 1] / kdata[i, j + 1]) ** (
- (T - Tlow) / (Thigh - Tlow))
+ klow = kdata[i, j] * (kdata[i + 1, j] / kdata[i, j]) ** ((T - Tlow) / (Thigh - Tlow)).real
+ khigh = kdata[i, j + 1] * (kdata[i + 1, j + 1] / kdata[i, j + 1]) ** ((T - Tlow) / (Thigh - Tlow)).real
k = klow * (khigh / klow) ** (log(P / Plow) / log(Phigh / Plow))
break
diff --git a/rmgpy/molecule/group.pxd b/rmgpy/molecule/group.pxd
index 83b4ae83c6..ada4072f78 100644
--- a/rmgpy/molecule/group.pxd
+++ b/rmgpy/molecule/group.pxd
@@ -31,6 +31,12 @@ cimport rmgpy.molecule.molecule as mol
from cpython cimport bool
################################################################################
+cdef bool check_set(super_list, sub_list)
+
+cdef list add_cb_atom_to_ring(ring, cb_atom)
+
+cdef list merge_overlapping_benzene_rings(ring1, ring2, od)
+
cdef class GroupAtom(Vertex):
cdef public list atomtype
diff --git a/rmgpy/molecule/group.py b/rmgpy/molecule/group.py
index 1dc56cc7f7..e9a6e6433f 100644
--- a/rmgpy/molecule/group.py
+++ b/rmgpy/molecule/group.py
@@ -46,8 +46,91 @@
from rmgpy.molecule.graph import Vertex, Edge, Graph
from rmgpy.molecule.fragment import CuttingLabel
+# helper functions
+# these were originall nested inside the indicated parent function, but when we upgraded to
+# Cython 3 this was no longer allowed - thus, they now live here.
-################################################################################
+# add_implicit_benzene
+def check_set(super_list, sub_list):
+ """
+ Args:
+ super_list: list to check if superset of partList
+ sub_list: list to check if subset of superList
+
+ Returns: Boolean to see if super_list is a superset of sub_list
+
+ """
+ super_set = set(super_list)
+ sub_set = set(sub_list)
+ return super_set.issuperset(sub_set)
+
+def add_cb_atom_to_ring(ring, cb_atom):
+ """
+ Every 'Cb' atom belongs in exactly one benzene ring. This function checks
+ adds the cb_atom to the ring (in connectivity order) if the cb_atom is connected
+ to any the last or first atom in the partial ring.
+
+ Args:
+ ring: list of :class:GroupAtoms representing a partial ring to merge
+ cb_atom: :class:GroupAtom with atomtype 'Cb'
+
+ Returns: If cb_atom connects to the beginning or end of ring, returns a
+ new list of the merged ring, otherwise an empty list
+
+ """
+
+ merged_ring = []
+ # ring already complete
+ if len(ring) == 6: return merged_ring
+ for atom2, bond12 in cb_atom.bonds.items():
+ if bond12.is_benzene():
+ if atom2 is ring[-1]:
+ merged_ring = ring + [cb_atom]
+ elif atom2 is ring[0]:
+ merged_ring = [cb_atom] + ring
+
+ return merged_ring
+
+def merge_overlapping_benzene_rings(ring1, ring2, od):
+ """
+ The input arguments of rings are always in the order that the atoms appear
+ inside the ring. That is, each atom is connected to the ones adjacent on the
+ list.
+
+ Args:
+ ring1: list of :class:GroupAtoms representing first partial ring to merge
+ ring2: list :class:GroupAtoms representing second partial ring to merge
+ od: in for overlap distance
+
+ This function tries to see if the beginning or ends of each list have the
+ same atom objects, i.e the two part rings should be merged together.
+
+ Returns: If rings are mergable, returns a new list of the merged ring, otherwise
+ an empty list
+
+ """
+ new_ring = []
+ # ring already complete
+ if len(ring1) == 6 or len(ring2) == 6: return new_ring
+
+ # start of ring1 matches end of ring2
+ match_list1 = [x1 is x2 for x1, x2 in zip(ring1[-od:], ring2[:od])]
+ # end of ring1 matches end of ring2
+ match_list2 = [x1 is x2 for x1, x2 in zip(ring1[-od:], ring2[:od - 1:-1])]
+ # start of ring1 matches end of ring2
+ match_list3 = [x1 is x2 for x1, x2 in zip(ring1[:od], ring2[-od:])]
+ # start of ring1 matches start of ring2
+ match_list4 = [x1 is x2 for x1, x2 in zip(ring1[:od], ring2[od::-1])]
+ if False not in match_list1:
+ new_ring = ring1 + ring2[od:]
+ elif False not in match_list2:
+ new_ring = ring1 + ring2[-od - 1::-1]
+ elif False not in match_list3:
+ new_ring = ring2[:-od] + ring1
+ elif False not in match_list4:
+ new_ring = ring2[:od - 1:-1] + ring1
+
+ return new_ring
class GroupAtom(Vertex):
"""
@@ -2531,89 +2614,6 @@ def add_implicit_benzene(self):
# Note that atomtypes like N5bd are mostly referred to as Cb in this code,
# which was first written for just carbon.
- # First define some helper functions
- def check_set(super_list, sub_list):
- """
- Args:
- super_list: list to check if superset of partList
- sub_list: list to check if subset of superList
-
- Returns: Boolean to see if super_list is a superset of sub_list
-
- """
- super_set = set(super_list)
- sub_set = set(sub_list)
- return super_set.issuperset(sub_set)
-
- def add_cb_atom_to_ring(ring, cb_atom):
- """
- Every 'Cb' atom belongs in exactly one benzene ring. This function checks
- adds the cb_atom to the ring (in connectivity order) if the cb_atom is connected
- to any the last or first atom in the partial ring.
-
- Args:
- ring: list of :class:GroupAtoms representing a partial ring to merge
- cb_atom: :class:GroupAtom with atomtype 'Cb'
-
- Returns: If cb_atom connects to the beginning or end of ring, returns a
- new list of the merged ring, otherwise an empty list
-
- """
-
- merged_ring = []
- # ring already complete
- if len(ring) == 6: return merged_ring
- for atom2, bond12 in cb_atom.bonds.items():
- if bond12.is_benzene():
- if atom2 is ring[-1]:
- merged_ring = ring + [cb_atom]
- elif atom2 is ring[0]:
- merged_ring = [cb_atom] + ring
-
- return merged_ring
-
- def merge_overlapping_benzene_rings(ring1, ring2, od):
- """
- The input arguments of rings are always in the order that the atoms appear
- inside the ring. That is, each atom is connected to the ones adjacent on the
- list.
-
- Args:
- ring1: list of :class:GroupAtoms representing first partial ring to merge
- ring2: list :class:GroupAtoms representing second partial ring to merge
- od: in for overlap distance
-
- This function tries to see if the beginning or ends of each list have the
- same atom objects, i.e the two part rings should be merged together.
-
- Returns: If rings are mergable, returns a new list of the merged ring, otherwise
- an empty list
-
- """
- new_ring = []
- # ring already complete
- if len(ring1) == 6 or len(ring2) == 6: return new_ring
-
- # start of ring1 matches end of ring2
- match_list1 = [x1 is x2 for x1, x2 in zip(ring1[-od:], ring2[:od])]
- # end of ring1 matches end of ring2
- match_list2 = [x1 is x2 for x1, x2 in zip(ring1[-od:], ring2[:od - 1:-1])]
- # start of ring1 matches end of ring2
- match_list3 = [x1 is x2 for x1, x2 in zip(ring1[:od], ring2[-od:])]
- # start of ring1 matches start of ring2
- match_list4 = [x1 is x2 for x1, x2 in zip(ring1[:od], ring2[od::-1])]
- if False not in match_list1:
- new_ring = ring1 + ring2[od:]
- elif False not in match_list2:
- new_ring = ring1 + ring2[-od - 1::-1]
- elif False not in match_list3:
- new_ring = ring2[:-od] + ring1
- elif False not in match_list4:
- new_ring = ring2[:od - 1:-1] + ring1
-
- return new_ring
-
- #######################################################################################
# start of main algorithm
copy_group = deepcopy(self)
"""
diff --git a/rmgpy/molecule/inchi.py b/rmgpy/molecule/inchi.py
index 59e49d8522..438b2b1532 100644
--- a/rmgpy/molecule/inchi.py
+++ b/rmgpy/molecule/inchi.py
@@ -30,6 +30,7 @@
import itertools
import re
import warnings
+from operator import itemgetter
import cython
from rdkit import Chem
@@ -614,7 +615,7 @@ def create_augmented_layers(mol):
atom_indices = [atom_indices.index(i + 1) for i, atom in enumerate(molcopy.atoms)]
# sort the atoms based on the order of the atom indices
- molcopy.atoms = [x for (y, x) in sorted(zip(atom_indices, molcopy.atoms), key=lambda pair: pair[0])]
+ molcopy.atoms = [x for (y, x) in sorted(zip(atom_indices, molcopy.atoms), key=itemgetter(0))]
ulayer = _create_u_layer(molcopy, auxinfo)
@@ -704,17 +705,15 @@ def _convert_3_atom_2_bond_path(start, mol):
"""
from rmgpy.data.kinetics.family import ReactionRecipe
- def is_valid(mol):
- """Check if total bond order of oxygen atoms is smaller than 4."""
-
- for at in mol.atoms:
- if at.number == 8:
- order = at.get_total_bond_order()
- not_correct = order >= 4
- if not_correct:
- return False
-
- return True
+ # Check if total bond order of oxygen atoms is smaller than 4
+ is_mol_valid = True
+ for at in mol.atoms:
+ if at.number == 8:
+ order = at.get_total_bond_order()
+ not_correct = order >= 4
+ if not_correct:
+ is_mol_valid = False
+ break
paths = pathfinder.find_allyl_end_with_charge(start)
@@ -739,7 +738,7 @@ def is_valid(mol):
end.charge += 1 if end.charge < 0 else -1
recipe.apply_forward(mol)
- if is_valid(mol):
+ if is_mol_valid:
# unlabel atoms so that they never cause trouble downstream
for i, at in enumerate(path[::2]):
at.label = ''
diff --git a/rmgpy/molecule/molecule.pxd b/rmgpy/molecule/molecule.pxd
index 335486f76f..d8bc24f4bb 100644
--- a/rmgpy/molecule/molecule.pxd
+++ b/rmgpy/molecule/molecule.pxd
@@ -36,6 +36,8 @@ from rmgpy.molecule.graph cimport Vertex, Edge, Graph
################################################################################
cdef dict bond_orders
+cdef tuple _skip_first(in_tuple)
+
cdef class Atom(Vertex):
cdef public Element element
diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py
index 3f647c336a..5ac1a48702 100644
--- a/rmgpy/molecule/molecule.py
+++ b/rmgpy/molecule/molecule.py
@@ -41,6 +41,7 @@
from collections import OrderedDict, defaultdict
from copy import deepcopy
from urllib.parse import quote
+from operator import attrgetter
import cython
import numpy as np
@@ -62,6 +63,10 @@
################################################################################
+# helper function for sorting
+def _skip_first(in_tuple):
+ return in_tuple[1:]
+
bond_orders = {'S': 1, 'D': 2, 'T': 3, 'B': 1.5}
globals().update({
@@ -2543,30 +2548,21 @@ def get_aromatic_rings(self, rings=None, save_order=False):
if rings is None:
rings = self.get_relevant_cycles()
- def filter_fused_rings(_rings):
- """
- Given a list of rings, remove ones which share more than 2 atoms.
- """
- cython.declare(toRemove=set, i=cython.int, j=cython.int, toRemoveSorted=list)
-
- if len(_rings) < 2:
- return _rings
-
+ # Remove rings that share more than 3 atoms, since they cannot be planar
+ cython.declare(toRemove=set, j=cython.int, toRemoveSorted=list)
+ if len(rings) < 2:
+ pass
+ else:
to_remove = set()
- for i, j in itertools.combinations(range(len(_rings)), 2):
- if len(set(_rings[i]) & set(_rings[j])) > 2:
+ for i, j in itertools.combinations(range(len(rings)), 2):
+ if len(set(rings[i]) & set(rings[j])) > 2:
to_remove.add(i)
to_remove.add(j)
to_remove_sorted = sorted(to_remove, reverse=True)
for i in to_remove_sorted:
- del _rings[i]
-
- return _rings
-
- # Remove rings that share more than 3 atoms, since they cannot be planar
- rings = filter_fused_rings(rings)
+ del rings[i]
# Only keep rings with exactly 6 atoms, since RMG can only handle aromatic benzene
rings = [ring for ring in rings if len(ring) == 6]
@@ -2702,7 +2698,7 @@ def get_deterministic_sssr(self):
tup = (vertex, get_vertex_connectivity_value(vertex), -origin_conn_dict[vertex])
root_candidates_tups.append(tup)
- root_vertex = sorted(root_candidates_tups, key=lambda tup0: tup0[1:], reverse=True)[0][0]
+ root_vertex = sorted(root_candidates_tups, key=_skip_first, reverse=True)[0][0]
# Get all cycles involving the root vertex
cycles = graph0.get_all_cycles(root_vertex)
@@ -2719,7 +2715,7 @@ def get_deterministic_sssr(self):
-sum([v.get_total_bond_order() for v in cycle0]))
cycle_candidate_tups.append(tup)
- cycle = sorted(cycle_candidate_tups, key=lambda tup0: tup0[1:])[0][0]
+ cycle = sorted(cycle_candidate_tups, key=_skip_first)[0][0]
cycle_list.append(cycle)
@@ -2803,8 +2799,8 @@ def is_identical(self, other, strict=True):
if atom_ids == other_ids:
# If the two molecules have the same indices, then they might be identical
# Sort the atoms by ID
- atom_list = sorted(self.atoms, key=lambda x: x.id)
- other_list = sorted(other.atoms, key=lambda x: x.id)
+ atom_list = sorted(self.atoms, key=attrgetter('id'))
+ other_list = sorted(other.atoms, key=attrgetter('id'))
# If matching atom indices gives a valid mapping, then the molecules are fully identical
mapping = {}
diff --git a/rmgpy/molecule/resonance.pxd b/rmgpy/molecule/resonance.pxd
index e3b01db841..a28c0aa1ad 100644
--- a/rmgpy/molecule/resonance.pxd
+++ b/rmgpy/molecule/resonance.pxd
@@ -25,9 +25,15 @@
# #
###############################################################################
+cimport numpy as cnp
+
from rmgpy.molecule.graph cimport Vertex, Edge, Graph
from rmgpy.molecule.molecule cimport Atom, Bond, Molecule
+cdef int _sum_atom_ids(atom_list)
+
+cdef tuple _tuplize_bond(bond)
+
cpdef list populate_resonance_algorithms(dict features=?)
cpdef dict analyze_molecule(Graph mol, bint save_order=?)
@@ -60,6 +66,4 @@ cpdef list generate_kekule_structure(Graph mol)
cpdef list generate_clar_structures(Graph mol, bint save_order=?)
-cpdef list _clar_optimization(Graph mol, list constraints=?, max_num=?, save_order=?)
-
-cpdef list _clar_transformation(Graph mol, list aromatic_ring)
+cpdef list _clar_optimization(Graph mol, list recursion_constraints=?, int clar_number=?, bint save_order=?)
diff --git a/rmgpy/molecule/resonance.py b/rmgpy/molecule/resonance.py
index b207568e71..3cb0bad9b8 100644
--- a/rmgpy/molecule/resonance.py
+++ b/rmgpy/molecule/resonance.py
@@ -52,17 +52,20 @@
"""
import logging
+from operator import attrgetter
import cython
+import numpy as np
+from scipy.optimize import Bounds, LinearConstraint, milp
import rmgpy.molecule.filtration as filtration
import rmgpy.molecule.pathfinder as pathfinder
-from rmgpy.exceptions import ILPSolutionError, KekulizationError, AtomTypeError, ResonanceError
+from rmgpy.exceptions import AtomTypeError, KekulizationError, ResonanceError
from rmgpy.molecule.adjlist import Saturator
+from rmgpy.molecule.fragment import CuttingLabel
from rmgpy.molecule.graph import Vertex
from rmgpy.molecule.kekulize import kekulize
from rmgpy.molecule.molecule import Atom, Bond, Molecule
-from rmgpy.molecule.fragment import CuttingLabel
def populate_resonance_algorithms(features=None):
@@ -905,8 +908,7 @@ def generate_clar_structures(mol, save_order=False):
Returns a list of :class:`Molecule` objects corresponding to the Clar structures.
"""
- cython.declare(output=list, mol_list=list, new_mol=Graph, aromatic_rings=list, bonds=list, solution=list,
- y=list, x=list, index=cython.int, bond=Edge, ring=list)
+ cython.declare(output=list, mol_list=list, new_mol=Graph, aromatic_rings=list, bonds=list, index=cython.int, bond=Edge, ring=list)
if not mol.is_cyclic():
return []
@@ -916,19 +918,18 @@ def generate_clar_structures(mol, save_order=False):
mol.assign_atom_ids()
try:
- output = _clar_optimization(mol, save_order=save_order)
- except ILPSolutionError:
+ solutions = _clar_optimization(mol, save_order=save_order)
+ except (RuntimeError, ValueError): # either a crash during optimization or the result was an empty tuple
# The optimization algorithm did not work on the first iteration
return []
mol_list = []
- for new_mol, aromatic_rings, bonds, solution in output:
-
+ for new_mol, aromatic_rings, bonds, solution in solutions:
# The solution includes a part corresponding to rings, y, and a part corresponding to bonds, x, using
# nomenclature from the paper. In y, 1 means the ring as a sextet, 0 means it does not.
# In x, 1 corresponds to a double bond, 0 either means a single bond or the bond is part of a sextet.
- y = solution[0:len(aromatic_rings)]
+ y = solution[:len(aromatic_rings)]
x = solution[len(aromatic_rings):]
# Apply results to molecule - double bond locations first
@@ -938,182 +939,202 @@ def generate_clar_structures(mol, save_order=False):
elif x[index] == 1:
bond.order = 2 # double
else:
- raise ValueError('Unaccepted bond value {0} obtained from optimization.'.format(x[index]))
+ raise ValueError(f'Unaccepted bond value {x[index]} obtained from optimization.')
# Then apply locations of aromatic sextets by converting to benzene bonds
for index, ring in enumerate(aromatic_rings):
if y[index] == 1:
- _clar_transformation(new_mol, ring)
+ for i, atom_1 in enumerate(ring):
+ for j, atom_2 in enumerate(ring):
+ if new_mol.has_bond(atom_1, atom_2):
+ new_mol.get_bond(atom_1, atom_2).order = 1.5
try:
new_mol.update_atomtypes()
except AtomTypeError:
pass
- else:
+ finally:
mol_list.append(new_mol)
return mol_list
+# helper functions for sorting
+def _sum_atom_ids(atom_list):
+ return sum(atom.id for atom in atom_list)
-def _clar_optimization(mol, constraints=None, max_num=None, save_order=False):
+def _tuplize_bond(bond):
+ return (bond.atom1.id, bond.atom2.id)
+
+
+def _clar_optimization(mol, recursion_constraints=None, clar_number=-1, save_order=False):
"""
- Implements linear programming algorithm for finding Clar structures. This algorithm maximizes the number
- of Clar sextets within the constraints of molecular geometry and atom valency.
+ Implements Mixed Integer Linear Programming for finding Clar structures.
+
+ First finds the Clar number (and an arbitrary structure with that number), then recursively
+ calls itself to enumerate more structures with that Clar number. No guarantees about
+ which structures will be found, or how many (we stop solving once solution would require
+ branching, which typically happens after at least a couple have already been found).
Returns a list of valid Clar solutions in the form of a tuple, with the following entries:
- [0] Molecule object
+ [0] Copy of mol
[1] List of aromatic rings
[2] List of bonds
- [3] Optimization solution
+ [3] Solution vector
- The optimization solution is a list of boolean values with sextet assignments followed by double bond assignments,
- with indices corresponding to the list of aromatic rings and list of bonds, respectively.
+ The solution vector is a binary integer list of length (# aromatic rings + # bonds) indicating
+ if each ring is aromatic and if each bond is double.
- Method adapted from:
- Hansen, P.; Zheng, M. The Clar Number of a Benzenoid Hydrocarbon and Linear Programming.
- J. Math. Chem. 1994, 15 (1), 93–107.
+ This implementation follows the original implementation very closely (Hansen, P.; Zheng, M. The
+ Clar Number of a Benzenoid Hydrocarbon and Linear Programming. J. Math. Chem. 1994, 15 (1), 93–107.)
+ with only slight modifications to prevent changing bond orders for exocyclic atoms.
"""
- cython.declare(molecule=Graph, aromatic_rings=list, exo=list, l=cython.int, m=cython.int, n=cython.int,
- a=list, objective=list, status=cython.int, solution=list, innerSolutions=list)
-
- from lpsolve55 import lpsolve
-
- # Make a copy of the molecule so we don't destroy the original
+ cython.declare(
+ molecule=Graph,
+ aromatic_rings=list,
+ exo_info=list,
+ n_ring=cython.int,
+ n_atom=cython.int,
+ n_bond=cython.int,
+ A=list,
+ solutions=list,
+ )
+
+ # after we find all the Clar structures we will modify the molecule bond orders to be aromatic,
+ # so we make explicit copies to avoid overwrites.
molecule = mol.copy(deep=True)
aromatic_rings = molecule.get_aromatic_rings(save_order=save_order)[0]
- aromatic_rings.sort(key=lambda x: sum([atom.id for atom in x]))
- if not aromatic_rings:
+ # doesn't work for system without aromatic rings, just exit
+ if len(aromatic_rings) == 0:
return []
- # Get list of atoms that are in rings
+ # stability between multiple runs
+ aromatic_rings.sort(key=_sum_atom_ids)
+
+ # Cython doesn't allow mutable defaults, so we just set it here instead
+ if recursion_constraints is None:
+ recursion_constraints = []
+
+ # Get set of atoms that are in rings
atoms = set()
for ring in aromatic_rings:
atoms.update(ring)
- atoms = sorted(atoms, key=lambda x: x.id)
+ atoms = sorted(atoms, key=attrgetter('id'))
- # Get list of bonds involving the ring atoms, ignoring bonds to hydrogen
+ # Get set of bonds involving the ring atoms, ignoring bonds to hydrogen
bonds = set()
for atom in atoms:
bonds.update([atom.bonds[key] for key in atom.bonds.keys() if key.is_non_hydrogen()])
- bonds = sorted(bonds, key=lambda x: (x.atom1.id, x.atom2.id))
+ bonds = sorted(bonds, key=_tuplize_bond)
- # Identify exocyclic bonds, and save their bond orders
- exo = []
+ # identify exocyclic bonds, and save their order if exocyclic
+ exo_info = []
for bond in bonds:
if bond.atom1 not in atoms or bond.atom2 not in atoms:
+ # save order for exocyclic
if bond.is_double():
- exo.append(1)
+ exo_info.append(1)
else:
- exo.append(0)
- else:
- exo.append(None)
+ exo_info.append(0)
+ else: # not exocyclic
+ exo_info.append(None)
# Dimensions
- l = len(aromatic_rings)
- m = len(atoms)
- n = l + len(bonds)
+ n_ring = len(aromatic_rings)
+ n_atom = len(atoms)
+ n_bond = len(bonds)
+
+ # The aromaticity assignment problem is formulated as an MILP problem
+ # minimize:
+ # c @ x
+ # such that
+ # b_l <= A @ x <= b_u
+ # l <= x <= u
+ # x are integers
# Connectivity matrix which indicates which rings and bonds each atom is in
- # Part of equality constraint Ax=b
- a = []
+ # Part of equality constraint A_eq @ x = b_eq
+ A = []
for atom in atoms:
in_ring = [1 if atom in ring else 0 for ring in aromatic_rings]
in_bond = [1 if atom in [bond.atom1, bond.atom2] else 0 for bond in bonds]
- a.append(in_ring + in_bond)
+ A.append(in_ring + in_bond)
+ initial_constraint = [LinearConstraint( # i.e. an equality constraint
+ A=np.array(A, dtype=int), lb=np.ones(n_atom, dtype=int), ub=np.ones(n_atom, dtype=int)
+ )]
+
+ # on recursive calls we already know the Clar number, so we can additionally constrain the system
+ # to only find structures with this number. Allows detecting when all formulae have been found and,
+ # in theory, should solve faster
+ if clar_number != -1:
+ initial_constraint += [
+ LinearConstraint(
+ A=np.array([1] * n_ring + [0] * n_bond),
+ ub=clar_number,
+ lb=clar_number,
+ ),
+ ]
# Objective vector for optimization: sextets have a weight of 1, double bonds have a weight of 0
- objective = [1] * l + [0] * len(bonds)
-
- # Solve LP problem using lpsolve
- lp = lpsolve('make_lp', m, n) # initialize lp with constraint matrix with m rows and n columns
- lpsolve('set_verbose', lp, 2) # reduce messages from lpsolve
- lpsolve('set_obj_fn', lp, objective) # set objective function
- lpsolve('set_maxim', lp) # set solver to maximize objective
- lpsolve('set_mat', lp, a) # set left hand side to constraint matrix
- lpsolve('set_rh_vec', lp, [1] * m) # set right hand side to 1 for all constraints
- for i in range(m): # set all constraints as equality constraints
- lpsolve('set_constr_type', lp, i + 1, '=')
- lpsolve('set_binary', lp, [True] * n) # set all variables to be binary
-
- # Constrain values of exocyclic bonds, since we don't want to modify them
- for i in range(l, n):
- if exo[i - l] is not None:
- # NOTE: lpsolve indexes from 1, so the variable we're changing should be i + 1
- lpsolve('set_bounds', lp, i + 1, exo[i - l], exo[i - l])
-
- # Add constraints to problem if provided
- if constraints is not None:
- for constraint in constraints:
- try:
- lpsolve('add_constraint', lp, constraint[0], '<=', constraint[1])
- except Exception as e:
- logging.debug('Unable to add constraint: {0} <= {1}'.format(constraint[0], constraint[1]))
- logging.debug(mol.to_adjacency_list())
- if str(e) == 'invalid vector.':
- raise ILPSolutionError('Unable to add constraint, likely due to '
- 'inconsistent aromatic ring perception.')
- else:
- raise
-
- status = lpsolve('solve', lp)
- obj_val, solution = lpsolve('get_solution', lp)[0:2]
- lpsolve('delete_lp', lp) # Delete the LP problem to clear up memory
+ # we negate because the original problem is formulated as maximization, whereas SciPy only does min
+ c = -np.array([1] * n_ring + [0] * n_bond, dtype=int)
+
+ # variable bounds
+ # - binary problem, so everything must be either zero or one
+ # - rings are simply 0 or 1
+ # - bonds are also 0 or 1, except exocyclic bonds which must remain unchanged
+ bounds = Bounds(
+ lb=np.array([0] * n_ring + [0 if b is None else b for b in exo_info], dtype=int),
+ ub=np.array([1] * n_ring + [1 if b is None else b for b in exo_info], dtype=int),
+ )
+
+ result = milp(
+ c=c,
+ integrality=1,
+ bounds=bounds,
+ constraints=initial_constraint + recursion_constraints,
+ options={'time_limit': 10},
+ )
+
+ if (status := result.status) != 0:
+ if status == 2: # infeasible
+ raise RuntimeError("All valid Clar formulae have been enumerated!")
+ else:
+ raise RuntimeError(f"Optimization failed (Exit Code {status}) for an unexpected reason '{result.message}'")
- # Check that optimization was successful
- if status != 0:
- raise ILPSolutionError('Optimization could not find a valid solution.')
+ _clar_number, solution = -result.fun, result.x
- # Check that we the result contains at least one aromatic sextet
- if obj_val == 0:
+ # optimization may have reached a bad local minimum - this case is rare
+ if _clar_number == 0:
return []
- # Check that the solution contains the maximum number of sextets possible
- if max_num is None:
- max_num = obj_val # This is the first solution, so the result should be an upper limit
- elif obj_val < max_num:
- raise ILPSolutionError('Optimization obtained a sub-optimal solution.')
-
+ # on later runs, non-integer solutions occur - branching might be able to find actual solutions,
+ # but we just call it good enough here
if any([x != 1 and x != 0 for x in solution]):
- raise ILPSolutionError('Optimization obtained a non-integer solution.')
-
- # Generate constraints based on the solution obtained
- y = solution[0:l]
- new_a = y + [0] * len(bonds)
- new_b = sum(y) - 1
- if constraints is not None:
- constraints.append((new_a, new_b))
- else:
- constraints = [(new_a, new_b)]
+ raise RuntimeError("Optimization obtained a non-integer solution - no more formulae will be enumerated.")
+
+ # first solution, so the result should be an upper limit
+ if clar_number == -1:
+ clar_number = _clar_number
+
+ selected_sextets = list(solution[:n_ring])
+ # restrict those rings which were selected for the current clar structure
+ # from forming it again by requiring that their sum is 1 less than it used
+ # to be
+ recursion_constraints += [
+ LinearConstraint(
+ A=np.array(selected_sextets + [0] * n_bond),
+ ub=clar_number - 1,
+ lb=0,
+ ),
+ ]
# Run optimization with additional constraints
try:
- inner_solutions = _clar_optimization(mol, constraints=constraints, max_num=max_num, save_order=save_order)
- except ILPSolutionError:
+ inner_solutions = _clar_optimization(mol, recursion_constraints=recursion_constraints, clar_number=clar_number, save_order=save_order)
+ except RuntimeError as e:
+ logging.debug(f"Clar Optimization stopped: {e}")
inner_solutions = []
return inner_solutions + [(molecule, aromatic_rings, bonds, solution)]
-
-
-def _clar_transformation(mol, aromatic_ring):
- """
- Performs Clar transformation for given ring in a molecule, ie. conversion to aromatic sextet.
-
- Args:
- mol a :class:`Molecule` object
- aromaticRing a list of :class:`Atom` objects corresponding to an aromatic ring in mol
-
- This function directly modifies the input molecule and does not return anything.
- """
- cython.declare(bondList=list, i=cython.int, atom1=Vertex, atom2=Vertex, bond=Edge)
-
- bond_list = []
-
- for i, atom1 in enumerate(aromatic_ring):
- for atom2 in aromatic_ring[i + 1:]:
- if mol.has_bond(atom1, atom2):
- bond_list.append(mol.get_bond(atom1, atom2))
-
- for bond in bond_list:
- bond.order = 1.5
diff --git a/rmgpy/molecule/util.py b/rmgpy/molecule/util.py
index d1093f4582..46f8c3504c 100644
--- a/rmgpy/molecule/util.py
+++ b/rmgpy/molecule/util.py
@@ -30,8 +30,8 @@
import itertools
import re
-from rmgpy.molecule.molecule import Molecule
from rmgpy.molecule.fragment import Fragment
+from rmgpy.molecule.molecule import Molecule
def get_element_count(obj):
@@ -47,7 +47,7 @@ def get_element_count(obj):
match = re.match(r"([a-z]+)([0-9]*)", piece, re.I)
if match:
element, count = match.groups()
- if count is '':
+ if count == '':
count = 1
if element in element_count:
element_count[element] += int(count)
diff --git a/rmgpy/pdep/reaction.pyx b/rmgpy/pdep/reaction.pyx
index 6a6327c5d7..d8e6942de4 100644
--- a/rmgpy/pdep/reaction.pyx
+++ b/rmgpy/pdep/reaction.pyx
@@ -308,7 +308,7 @@ def apply_inverse_laplace_transform_method(transition_state,
and abs(dens_states[r - m, s]) > 1e-12 \
and abs(dens_states[r - m - 1, s]) > 1e-12:
num = dens_states[r - m, s] * (dens_states[r - m - 1, s] / dens_states[r - m, s]) \
- ** (-rem / (e_list[r - m - 1] - e_list[r - m]))
+ ** (-rem / (e_list[r - m - 1] - e_list[r - m])).real
k[r, s] = freq_factor * num / dens_states[r, s]
elif n >= n_crit:
diff --git a/rmgpy/pdep/sls.py b/rmgpy/pdep/sls.py
index b183a0bba8..12a0dd8e05 100644
--- a/rmgpy/pdep/sls.py
+++ b/rmgpy/pdep/sls.py
@@ -31,20 +31,31 @@
Contains functionality for directly simulating the master equation
and implementing the SLS master equation reduction method
"""
+
+import logging
import sys
-from diffeqpy import de
-from julia import Main
-import scipy.sparse as sparse
import numpy as np
import scipy.linalg
-import mpmath
import scipy.optimize as opt
+import scipy.sparse as sparse
import rmgpy.constants as constants
-from rmgpy.pdep.reaction import calculate_microcanonical_rate_coefficient
from rmgpy.pdep.me import generate_full_me_matrix, states_to_configurations
-from rmgpy.statmech.translation import IdealGasTranslation
+from rmgpy.rmg.reactionmechanismsimulator_reactors import to_julia
+
+NO_JULIA = False
+try:
+ from juliacall import Main
+
+ Main.seval("using ReactionMechanismSimulator.SciMLBase")
+ Main.seval("using ReactionMechanismSimulator.Sundials")
+except Exception as e:
+ logging.info(
+ f"Unable to import Julia dependencies, original error: {str(e)}"
+ ". Master equation method 'ode' will not be available on this execution."
+ )
+ NO_JULIA = True
def get_initial_condition(network, x0, indices):
@@ -85,30 +96,36 @@ def get_initial_condition(network, x0, indices):
def solve_me(M, p0, t):
- f = Main.eval(
+ f = Main.seval(
"""
-function f(u, M, t)
- return M*u
+function f(du, u, M, t)
+ du .= M * u
+ return du
end"""
)
- jac = Main.eval(
+ jac = Main.seval(
"""
-function jac(u, M, t)
- return M
+function jac(J, u, M, t)
+ J .= M
+ return J
end"""
)
+ p0 = to_julia(p0)
+ M = to_julia(M)
tspan = (0.0, t)
- fcn = de.ODEFunction(f, jac=jac)
- prob = de.ODEProblem(fcn, p0, tspan, M)
- sol = de.solve(prob, solver=de.CVODE_BDF(), abstol=1e-16, reltol=1e-6)
+ fcn = Main.ODEFunction(f, jac=jac)
+ prob = Main.ODEProblem(fcn, p0, tspan, M)
+ sol = Main.solve(prob, Main.CVODE_BDF(), abstol=1e-16, reltol=1e-6)
return sol
def solve_me_fcns(f, jac, M, p0, t):
+ p0 = to_julia(p0)
+ M = to_julia(M)
tspan = (0.0, t)
- fcn = de.ODEFunction(f, jac=jac)
- prob = de.ODEProblem(fcn, p0, tspan, M)
- sol = de.solve(prob, solver=de.CVODE_BDF(), abstol=1e-16, reltol=1e-6)
+ fcn = Main.ODEFunction(f, jac=jac)
+ prob = Main.ODEProblem(fcn, p0, tspan, M)
+ sol = Main.solve(prob, Main.CVODE_BDF(), abstol=1e-16, reltol=1e-6)
return sol
@@ -150,16 +167,23 @@ def get_rate_coefficients_SLS(network, T, P, method="mexp", neglect_high_energy_
tau = np.abs(1.0 / fastest_reaction)
if method == "ode":
- f = Main.eval(
+ if NO_JULIA:
+ raise RuntimeError(
+ "Required Julia dependencies for method 'ode' are not installed.\n"
+ "Please follow the steps to install Julia dependencies at https://reactionmechanismgenerator.github.io/RMG-Py/users/rmg/installation/anacondaDeveloper.html."
+ )
+ f = Main.seval(
"""
-function f(u,M,t)
- return M*u
+function f(du, u, M, t)
+ du .= M * u
+ return du
end"""
)
- jac = Main.eval(
+ jac = Main.seval(
"""
-function jac(u,M,t)
- return M
+function jac(J, u, M, t)
+ J .= M
+ return J
end"""
)
@@ -247,7 +271,7 @@ def run_equilibrium(network, channel):
xseq = []
dxdtseq = []
- # Single domainant source simulations
+ # Single dominant source simulations
for i, isomer in enumerate(isomers):
xsout, dxdtout = run_single_source(network, isomer)
for j in range(len(xsout)):
diff --git a/rmgpy/reaction.pxd b/rmgpy/reaction.pxd
index 15feecef3f..c6c8b35601 100644
--- a/rmgpy/reaction.pxd
+++ b/rmgpy/reaction.pxd
@@ -38,6 +38,8 @@ cimport numpy as np
################################################################################
+cdef tuple get_sorting_key(spc)
+
cdef class Reaction:
cdef public int index
cdef public str label
diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py
index b8ab3951c8..1997824afc 100644
--- a/rmgpy/reaction.py
+++ b/rmgpy/reaction.py
@@ -43,7 +43,7 @@
import math
import os.path
from copy import deepcopy
-from functools import reduce
+from functools import reduce, partial
from urllib.parse import quote
import cython
@@ -65,6 +65,16 @@
################################################################################
+# helper function for sorting
+def get_sorting_key(spc):
+ # List of elements to sort by, order is intentional
+ numbers = [6, 8, 7, 14, 16, 15, 17, 53, 9, 35] # C, O, N, Si, S, P, Cl, I, F, Br
+ ele_count = dict([(n,0) for n in numbers])
+ for atom in spc.molecule[0].atoms:
+ if isinstance(atom, Atom) and atom.element.number in numbers:
+ ele_count[atom.element.number] += 1
+ return tuple(ele_count[n] for n in numbers)
+
class Reaction:
"""
@@ -1516,14 +1526,7 @@ def generate_pairs(self):
reactants = self.reactants[:]
products = self.products[:]
- def get_sorting_key(spc):
- # List of elements to sort by, order is intentional
- numbers = [6, 8, 7, 14, 16, 15, 17, 53, 9, 35] # C, O, N, Si, S, P, Cl, I, F, Br
- ele_count = dict([(n,0) for n in numbers])
- for atom in spc.molecule[0].atoms:
- if isinstance(atom, Atom) and atom.element.number in numbers:
- ele_count[atom.element.number] += 1
- return tuple(ele_count[n] for n in numbers)
+
# Sort the reactants and products by element counts
reactants.sort(key=get_sorting_key)
products.sort(key=get_sorting_key)
@@ -1769,7 +1772,7 @@ def get_reduced_mass(self, reverse=False):
mass_list = [spc.molecule[0].get_molecular_weight() for spc in self.products]
else:
mass_list = [spc.molecule[0].get_molecular_weight() for spc in self.reactants]
- reduced_mass = reduce((lambda x, y: x * y), mass_list) / sum(mass_list)
+ reduced_mass = np.prod(mass_list) / sum(mass_list)
return reduced_mass
def get_mean_sigma_and_epsilon(self, reverse=False):
@@ -1795,7 +1798,8 @@ def get_mean_sigma_and_epsilon(self, reverse=False):
if any([x == 0 for x in sigmas + epsilons]):
raise ValueError
mean_sigmas = sum(sigmas) / num_of_spcs
- mean_epsilons = reduce((lambda x, y: x * y), epsilons) ** (1 / len(epsilons))
+ mean_epsilons = np.prod(epsilons) ** (1 / len(epsilons))
+
return mean_sigmas, mean_epsilons
def generate_high_p_limit_kinetics(self):
@@ -1805,6 +1809,15 @@ def generate_high_p_limit_kinetics(self):
"""
raise NotImplementedError("generate_high_p_limit_kinetics is not implemented for all Reaction subclasses.")
+def _same_object(object1, object2, _check_identical=False, _only_check_label=False,
+ _generate_initial_map=False, _strict=True, _save_order=False):
+ if _only_check_label:
+ return str(object1) == str(object2)
+ elif _check_identical:
+ return object1.is_identical(object2, strict=_strict)
+ else:
+ return object1.is_isomorphic(object2, generate_initial_map=_generate_initial_map,
+ strict=_strict, save_order=_save_order)
def same_species_lists(list1, list2, check_identical=False, only_check_label=False, generate_initial_map=False,
strict=True, save_order=False):
@@ -1825,46 +1838,45 @@ def same_species_lists(list1, list2, check_identical=False, only_check_label=Fal
Returns:
``True`` if the lists are the same and ``False`` otherwise
"""
-
- def same(object1, object2, _check_identical=check_identical, _only_check_label=only_check_label,
- _generate_initial_map=generate_initial_map, _strict=strict, _save_order=save_order):
- if _only_check_label:
- return str(object1) == str(object2)
- elif _check_identical:
- return object1.is_identical(object2, strict=_strict)
- else:
- return object1.is_isomorphic(object2, generate_initial_map=_generate_initial_map,
- strict=_strict, save_order=_save_order)
+
+ same_object_passthrough = partial(
+ _same_object,
+ _check_identical=check_identical,
+ _only_check_label=only_check_label,
+ _generate_initial_map=generate_initial_map,
+ _strict=strict,
+ _save_order=save_order,
+ )
if len(list1) == len(list2) == 1:
- if same(list1[0], list2[0]):
+ if same_object_passthrough(list1[0], list2[0]):
return True
elif len(list1) == len(list2) == 2:
- if same(list1[0], list2[0]) and same(list1[1], list2[1]):
+ if same_object_passthrough(list1[0], list2[0]) and same_object_passthrough(list1[1], list2[1]):
return True
- elif same(list1[0], list2[1]) and same(list1[1], list2[0]):
+ elif same_object_passthrough(list1[0], list2[1]) and same_object_passthrough(list1[1], list2[0]):
return True
elif len(list1) == len(list2) == 3:
- if same(list1[0], list2[0]):
- if same(list1[1], list2[1]):
- if same(list1[2], list2[2]):
+ if same_object_passthrough(list1[0], list2[0]):
+ if same_object_passthrough(list1[1], list2[1]):
+ if same_object_passthrough(list1[2], list2[2]):
return True
- elif same(list1[1], list2[2]):
- if same(list1[2], list2[1]):
+ elif same_object_passthrough(list1[1], list2[2]):
+ if same_object_passthrough(list1[2], list2[1]):
return True
- elif same(list1[0], list2[1]):
- if same(list1[1], list2[0]):
- if same(list1[2], list2[2]):
+ elif same_object_passthrough(list1[0], list2[1]):
+ if same_object_passthrough(list1[1], list2[0]):
+ if same_object_passthrough(list1[2], list2[2]):
return True
- elif same(list1[1], list2[2]):
- if same(list1[2], list2[0]):
+ elif same_object_passthrough(list1[1], list2[2]):
+ if same_object_passthrough(list1[2], list2[0]):
return True
- elif same(list1[0], list2[2]):
- if same(list1[1], list2[0]):
- if same(list1[2], list2[1]):
+ elif same_object_passthrough(list1[0], list2[2]):
+ if same_object_passthrough(list1[1], list2[0]):
+ if same_object_passthrough(list1[2], list2[1]):
return True
- elif same(list1[1], list2[1]):
- if same(list1[2], list2[0]):
+ elif same_object_passthrough(list1[1], list2[1]):
+ if same_object_passthrough(list1[2], list2[0]):
return True
elif len(list1) == len(list2):
raise NotImplementedError("Can't check isomorphism of lists with {0} species/molecules".format(len(list1)))
diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py
index 117c19644c..ba48ceb48a 100644
--- a/rmgpy/rmg/input.py
+++ b/rmgpy/rmg/input.py
@@ -35,23 +35,36 @@
from rmgpy import settings
from rmgpy.data.base import Entry
+from rmgpy.data.solvation import SolventData
+from rmgpy.data.surface import MetalDatabase
+from rmgpy.data.vaporLiquidMassTransfer import (
+ liquidVolumetricMassTransferCoefficientPowerLaw,
+)
from rmgpy.exceptions import DatabaseError, InputError
from rmgpy.molecule import Molecule
+from rmgpy.molecule.fragment import Fragment
from rmgpy.molecule.group import Group
-from rmgpy.quantity import Quantity, Energy, RateCoefficient, SurfaceConcentration
+from rmgpy.quantity import Energy, Quantity, RateCoefficient, SurfaceConcentration
from rmgpy.rmg.model import CoreEdgeReactionModel
+from rmgpy.rmg.reactionmechanismsimulator_reactors import (
+ NO_JULIA,
+ ConstantTLiquidSurfaceReactor,
+ ConstantTPIdealGasReactor,
+ ConstantTVLiquidReactor,
+ ConstantVIdealGasReactor,
+ Reactor,
+)
from rmgpy.rmg.settings import ModelSettings, SimulatorSettings
-from rmgpy.solver.termination import TerminationTime, TerminationConversion, TerminationRateRatio
from rmgpy.solver.liquid import LiquidReactor
from rmgpy.solver.mbSampled import MBSampledReactor
from rmgpy.solver.simple import SimpleReactor
from rmgpy.solver.surface import SurfaceReactor
+from rmgpy.solver.termination import (
+ TerminationConversion,
+ TerminationRateRatio,
+ TerminationTime,
+)
from rmgpy.util import as_list
-from rmgpy.data.surface import MetalDatabase
-from rmgpy.rmg.reactors import Reactor, ConstantVIdealGasReactor, ConstantTLiquidSurfaceReactor, ConstantTVLiquidReactor, ConstantTPIdealGasReactor
-from rmgpy.data.vaporLiquidMassTransfer import liquidVolumetricMassTransferCoefficientPowerLaw
-from rmgpy.molecule.fragment import Fragment
-from rmgpy.data.solvation import SolventData
################################################################################
@@ -1580,6 +1593,11 @@ def read_input_file(path, rmg0):
exec(f.read(), global_context, local_context)
except (NameError, TypeError, SyntaxError) as e:
logging.error('The input file "{0}" was invalid:'.format(full_path))
+ if NO_JULIA:
+ logging.error(
+ "During runtime, import of Julia dependencies failed. To use phase systems and RMS reactors, install RMG-Py with RMS."
+ " (https://reactionmechanismgenerator.github.io/RMG-Py/users/rmg/installation/anacondaDeveloper.html)"
+ )
logging.exception(e)
raise
finally:
diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py
index fb4ae33ef8..3f4b574c77 100644
--- a/rmgpy/rmg/main.py
+++ b/rmgpy/rmg/main.py
@@ -52,8 +52,6 @@
from scipy.optimize import brute
import rmgpy.util as util
-from rmgpy.rmg.model import Species, CoreEdgeReactionModel
-from rmgpy.rmg.pdep import PDepNetwork
from rmgpy import settings
from rmgpy.chemkin import ChemkinWriter
from rmgpy.constraints import fails_species_constraints
@@ -61,26 +59,32 @@
from rmgpy.data.kinetics.family import TemplateReaction
from rmgpy.data.kinetics.library import KineticsLibrary, LibraryReaction
from rmgpy.data.rmg import RMGDatabase
-from rmgpy.exceptions import ForbiddenStructureException, DatabaseError, CoreError, InputError
-from rmgpy.kinetics.diffusionLimited import diffusion_limiter
from rmgpy.data.vaporLiquidMassTransfer import vapor_liquid_mass_transfer
-from rmgpy.kinetics import ThirdBody
-from rmgpy.kinetics import Troe
+from rmgpy.exceptions import (
+ CoreError,
+ DatabaseError,
+ ForbiddenStructureException,
+ InputError,
+)
+from rmgpy.kinetics import ThirdBody, Troe
+from rmgpy.kinetics.diffusionLimited import diffusion_limiter
from rmgpy.molecule import Molecule
from rmgpy.qm.main import QMDatabaseWriter
from rmgpy.reaction import Reaction
-from rmgpy.rmg.listener import SimulationProfileWriter, SimulationProfilePlotter
+from rmgpy.rmg.listener import SimulationProfilePlotter, SimulationProfileWriter
+from rmgpy.rmg.model import CoreEdgeReactionModel, Species
from rmgpy.rmg.output import OutputHTMLWriter
-from rmgpy.rmg.pdep import PDepReaction
+from rmgpy.rmg.pdep import PDepNetwork, PDepReaction
+from rmgpy.rmg.reactionmechanismsimulator_reactors import NO_JULIA
+from rmgpy.rmg.reactionmechanismsimulator_reactors import Reactor as RMSReactor
from rmgpy.rmg.settings import ModelSettings
-from rmgpy.solver.base import TerminationTime, TerminationConversion
+from rmgpy.solver.base import TerminationConversion, TerminationTime
from rmgpy.solver.simple import SimpleReactor
from rmgpy.stats import ExecutionStatsWriter
from rmgpy.thermo.thermoengine import submit
from rmgpy.tools.plot import plot_sensitivity
from rmgpy.tools.uncertainty import Uncertainty, process_local_results
from rmgpy.yml import RMSWriter
-from rmgpy.rmg.reactors import Reactor
################################################################################
@@ -271,12 +275,6 @@ def load_input(self, path=None):
if self.solvent:
self.reaction_model.solvent_name = self.solvent
- if self.surface_site_density:
- self.reaction_model.surface_site_density = self.surface_site_density
- self.reaction_model.core.phase_system.phases["Surface"].site_density = self.surface_site_density.value_si
- self.reaction_model.edge.phase_system.phases["Surface"].site_density = self.surface_site_density.value_si
- self.reaction_model.coverage_dependence = self.coverage_dependence
-
self.reaction_model.verbose_comments = self.verbose_comments
self.reaction_model.save_edge_species = self.save_edge_species
@@ -513,6 +511,19 @@ def initialize(self, **kwargs):
# Read input file
self.load_input(self.input_file)
+
+ # Check if ReactionMechanismSimulator reactors are being used
+ # if RMS is not installed but the user attempted to use it, the load_input_file would have failed
+ # if RMS is installed but they did not use it, we can avoid extra work
+ # if RMS is not installed and they did not use it, we avoid calling certain functions that would raise an error
+ requires_rms = any(isinstance(reactor_system, RMSReactor) for reactor_system in self.reaction_systems)
+
+ if self.surface_site_density:
+ self.reaction_model.surface_site_density = self.surface_site_density
+ if requires_rms:
+ self.reaction_model.core.phase_system.phases["Surface"].site_density = self.surface_site_density.value_si
+ self.reaction_model.edge.phase_system.phases["Surface"].site_density = self.surface_site_density.value_si
+ self.reaction_model.coverage_dependence = self.coverage_dependence
if kwargs.get("restart", ""):
import rmgpy.rmg.input
@@ -556,7 +567,7 @@ def initialize(self, **kwargs):
self.load_database()
for reaction_system in self.reaction_systems:
- if isinstance(reaction_system, Reactor):
+ if isinstance(reaction_system, RMSReactor):
reaction_system.finish_termination_criteria()
# Load restart seed mechanism (if specified)
@@ -604,8 +615,9 @@ def initialize(self, **kwargs):
longDesc='',
)
solvent_mol = False
- self.reaction_model.core.phase_system.phases["Default"].set_solvent(solvent_data)
- self.reaction_model.edge.phase_system.phases["Default"].set_solvent(solvent_data)
+ if requires_rms:
+ self.reaction_model.core.phase_system.phases["Default"].set_solvent(solvent_data)
+ self.reaction_model.edge.phase_system.phases["Default"].set_solvent(solvent_data)
diffusion_limiter.enable(solvent_data, self.database.solvation)
logging.info("Setting solvent data for {0}".format(self.solvent))
@@ -647,12 +659,12 @@ def initialize(self, **kwargs):
# Seed mechanisms: add species and reactions from seed mechanism
# DON'T generate any more reactions for the seed species at this time
for seed_mechanism in self.seed_mechanisms:
- self.reaction_model.add_seed_mechanism_to_core(seed_mechanism, react=False)
+ self.reaction_model.add_seed_mechanism_to_core(seed_mechanism, react=False, requires_rms=requires_rms)
# Reaction libraries: add species and reactions from reaction library to the edge so
# that RMG can find them if their rates are large enough
for library, option in self.reaction_libraries:
- self.reaction_model.add_reaction_library_to_edge(library)
+ self.reaction_model.add_reaction_library_to_edge(library, requires_rms=requires_rms)
# Also always add in a few bath gases (since RMG-Java does)
for label, smiles in [("Ar", "[Ar]"), ("He", "[He]"), ("Ne", "[Ne]"), ("N2", "N#N")]:
@@ -728,35 +740,37 @@ def initialize(self, **kwargs):
# This is necessary so that the PDep algorithm can identify the bath gas
for spec in self.initial_species:
if not spec.reactive:
- self.reaction_model.enlarge(spec)
+ self.reaction_model.enlarge(spec, requires_rms=requires_rms)
for spec in self.initial_species:
if spec.reactive:
- self.reaction_model.enlarge(spec)
+ self.reaction_model.enlarge(spec, requires_rms=requires_rms)
# chatelak: store constant SPC indices in the reactor attributes if any constant SPC provided in the input file
# advantages to write it here: this is run only once (as species indexes does not change over the generation)
if self.solvent is not None:
for index, reaction_system in enumerate(self.reaction_systems):
if (
- not isinstance(reaction_system, Reactor) and reaction_system.const_spc_names is not None
- ): # if no constant species provided do nothing
+ not isinstance(reaction_system, RMSReactor)
+ ) and reaction_system.const_spc_names is not None: # if no constant species provided do nothing
reaction_system.get_const_spc_indices(self.reaction_model.core.species) # call the function to identify indices in the solver
self.initialize_reaction_threshold_and_react_flags()
if self.filter_reactions and self.init_react_tuples:
- self.react_init_tuples()
+ self.react_init_tuples(requires_rms=requires_rms)
self.reaction_model.initialize_index_species_dict()
self.initialize_seed_mech()
+ return requires_rms
- def register_listeners(self):
+ def register_listeners(self, requires_rms=False):
"""
Attaches listener classes depending on the options
found in the RMG input file.
"""
self.attach(ChemkinWriter(self.output_directory))
- self.attach(RMSWriter(self.output_directory))
+ if requires_rms:
+ self.attach(RMSWriter(self.output_directory))
if self.generate_output_html:
self.attach(OutputHTMLWriter(self.output_directory))
@@ -768,7 +782,7 @@ def register_listeners(self):
if self.save_simulation_profiles:
for index, reaction_system in enumerate(self.reaction_systems):
- if isinstance(reaction_system, Reactor):
+ if requires_rms and isinstance(reaction_system, RMSReactor):
typ = type(reaction_system)
raise InputError(f"save_simulation_profiles=True not compatible with reactor of type {typ}")
reaction_system.attach(SimulationProfileWriter(self.output_directory, index, self.reaction_model.core.species))
@@ -782,10 +796,10 @@ def execute(self, initialize=True, **kwargs):
"""
if initialize:
- self.initialize(**kwargs)
+ requires_rms = self.initialize(**kwargs)
# register listeners
- self.register_listeners()
+ self.register_listeners(requires_rms=requires_rms)
self.done = False
@@ -812,7 +826,7 @@ def execute(self, initialize=True, **kwargs):
# Update react flags
if self.filter_reactions:
# Run the reaction system to update threshold and react flags
- if isinstance(reaction_system, Reactor):
+ if requires_rms and isinstance(reaction_system, RMSReactor):
self.update_reaction_threshold_and_react_flags(
rxn_sys_unimol_threshold=np.zeros((len(self.reaction_model.core.species),), bool),
rxn_sys_bimol_threshold=np.zeros((len(self.reaction_model.core.species), len(self.reaction_model.core.species)), bool),
@@ -855,6 +869,7 @@ def execute(self, initialize=True, **kwargs):
unimolecular_react=self.unimolecular_react,
bimolecular_react=self.bimolecular_react,
trimolecular_react=self.trimolecular_react,
+ requires_rms=requires_rms,
)
if not np.isinf(self.model_settings_list[0].thermo_tol_keep_spc_in_edge):
@@ -867,7 +882,7 @@ def execute(self, initialize=True, **kwargs):
)
if not np.isinf(self.model_settings_list[0].thermo_tol_keep_spc_in_edge):
- self.reaction_model.thermo_filter_down(maximum_edge_species=self.model_settings_list[0].maximum_edge_species)
+ self.reaction_model.thermo_filter_down(maximum_edge_species=self.model_settings_list[0].maximum_edge_species, requires_rms=requires_rms)
logging.info("Completed initial enlarge edge step.\n")
@@ -933,7 +948,7 @@ def execute(self, initialize=True, **kwargs):
prune = False
try:
- if isinstance(reaction_system, Reactor):
+ if requires_rms and isinstance(reaction_system, RMSReactor):
(
terminated,
resurrected,
@@ -1026,7 +1041,7 @@ def execute(self, initialize=True, **kwargs):
# Add objects to enlarge to the core first
for objectToEnlarge in objects_to_enlarge:
- self.reaction_model.enlarge(objectToEnlarge)
+ self.reaction_model.enlarge(objectToEnlarge, requires_rms=requires_rms)
if model_settings.filter_reactions:
# Run a raw simulation to get updated reaction system threshold values
@@ -1035,7 +1050,7 @@ def execute(self, initialize=True, **kwargs):
temp_model_settings.tol_keep_in_edge = 0
if not resurrected:
try:
- if isinstance(reaction_system, Reactor):
+ if requires_rms and isinstance(reaction_system, RMSReactor):
(
terminated,
resurrected,
@@ -1104,7 +1119,7 @@ def execute(self, initialize=True, **kwargs):
skip_update=True,
)
logging.warning(
- "Reaction thresholds/flags for Reaction System {0} was not updated due " "to resurrection".format(index + 1)
+ "Reaction thresholds/flags for Reaction System {0} was not updated due to resurrection".format(index + 1)
)
logging.info("")
@@ -1127,13 +1142,14 @@ def execute(self, initialize=True, **kwargs):
unimolecular_react=self.unimolecular_react,
bimolecular_react=self.bimolecular_react,
trimolecular_react=self.trimolecular_react,
+ requires_rms=requires_rms,
)
if old_edge_size != len(self.reaction_model.edge.reactions) or old_core_size != len(self.reaction_model.core.reactions):
reactor_done = False
if not np.isinf(self.model_settings_list[0].thermo_tol_keep_spc_in_edge):
- self.reaction_model.thermo_filter_down(maximum_edge_species=model_settings.maximum_edge_species)
+ self.reaction_model.thermo_filter_down(maximum_edge_species=model_settings.maximum_edge_species, requires_rms=requires_rms)
max_num_spcs_hit = len(self.reaction_model.core.species) >= model_settings.max_num_species
@@ -1160,6 +1176,7 @@ def execute(self, initialize=True, **kwargs):
model_settings.tol_move_to_core,
model_settings.maximum_edge_species,
model_settings.min_species_exist_iterations_for_prune,
+ requires_rms=requires_rms,
)
# Perform garbage collection after pruning
collected = gc.collect()
@@ -1274,12 +1291,16 @@ def run_uncertainty_analysis(self):
if self.uncertainty is not None and self.uncertainty["global"]:
try:
import muq
- except ImportError:
- logging.error("Unable to import MUQ. Skipping global uncertainty analysis.")
+ except ImportError as ie:
+ logging.error(
+ f"Skipping global uncertainty analysis! Unable to import MUQ (original error: {str(ie)})."
+ "Install muq with 'conda install -c conda-forge muq'."
+ )
self.uncertainty["global"] = False
else:
- import re
import random
+ import re
+
from rmgpy.tools.canteramodel import Cantera
from rmgpy.tools.globaluncertainty import ReactorPCEFactory
@@ -1924,7 +1945,7 @@ def initialize_reaction_threshold_and_react_flags(self):
if self.trimolecular:
self.trimolecular_react[:num_restart_spcs, :num_restart_spcs, :num_restart_spcs] = False
- def react_init_tuples(self):
+ def react_init_tuples(self, requires_rms=False):
"""
Reacts tuples given in the react block
"""
@@ -1957,6 +1978,7 @@ def react_init_tuples(self):
unimolecular_react=self.unimolecular_react,
bimolecular_react=self.bimolecular_react,
trimolecular_react=self.trimolecular_react,
+ requires_rms=requires_rms,
)
def update_reaction_threshold_and_react_flags(
@@ -2251,7 +2273,7 @@ def __init__(self, reaction_system, bspc):
if isinstance(value, list):
self.Ranges[key] = [v.value_si for v in value]
- if isinstance(reaction_system, Reactor):
+ if isinstance(reaction_system, RMSReactor):
self.tmax = reaction_system.tf
else:
for term in reaction_system.termination:
@@ -2469,7 +2491,16 @@ def make_profile_graph(stats_file, force_graph_generation=False):
if display_found or force_graph_generation:
try:
- from gprof2dot import PstatsParser, DotWriter, SAMPLES, themes, TIME, TIME_RATIO, TOTAL_TIME, TOTAL_TIME_RATIO
+ from gprof2dot import (
+ SAMPLES,
+ TIME,
+ TIME_RATIO,
+ TOTAL_TIME,
+ TOTAL_TIME_RATIO,
+ DotWriter,
+ PstatsParser,
+ themes,
+ )
except ImportError:
logging.warning("Trouble importing from package gprof2dot. Unable to create a graph of the profile statistics.")
logging.warning("Try getting the latest version with something like `pip install --upgrade gprof2dot`.")
diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py
index ee6795806b..713e7b6ea0 100644
--- a/rmgpy/rmg/model.py
+++ b/rmgpy/rmg/model.py
@@ -44,21 +44,27 @@
from rmgpy.data.kinetics.depository import DepositoryReaction
from rmgpy.data.kinetics.family import KineticsFamily, TemplateReaction
from rmgpy.data.kinetics.library import KineticsLibrary, LibraryReaction
-from rmgpy.data.vaporLiquidMassTransfer import vapor_liquid_mass_transfer
-from rmgpy.molecule.group import Group
from rmgpy.data.rmg import get_db
+from rmgpy.data.vaporLiquidMassTransfer import vapor_liquid_mass_transfer
from rmgpy.display import display
from rmgpy.exceptions import ForbiddenStructureException
-from rmgpy.kinetics import KineticsData, Arrhenius
+from rmgpy.kinetics import Arrhenius, KineticsData
+from rmgpy.molecule.fragment import Fragment
+from rmgpy.molecule.group import Group
from rmgpy.quantity import Quantity
from rmgpy.reaction import Reaction
-from rmgpy.rmg.pdep import PDepReaction, PDepNetwork
+from rmgpy.rmg.decay import decay_species
+from rmgpy.rmg.pdep import PDepNetwork, PDepReaction
from rmgpy.rmg.react import react_all
+from rmgpy.rmg.reactionmechanismsimulator_reactors import (
+ NO_JULIA,
+ Interface,
+ Phase,
+ PhaseSystem,
+)
+from rmgpy.rmg.reactionmechanismsimulator_reactors import Reactor as RMSReactor
from rmgpy.species import Species
from rmgpy.thermo.thermoengine import submit
-from rmgpy.rmg.decay import decay_species
-from rmgpy.rmg.reactors import PhaseSystem, Phase, Interface, Reactor
-from rmgpy.molecule.fragment import Fragment
################################################################################
@@ -75,7 +81,7 @@ def __init__(self, species=None, reactions=None, phases=None, interfaces={}):
if phases is None:
phases = {"Default": Phase(), "Surface": Phase()}
interfaces = {frozenset({"Default", "Surface"}): Interface(list(phases.values()))}
- self.phase_system = PhaseSystem(phases, interfaces)
+ self.phase_system = None if NO_JULIA else PhaseSystem(phases, interfaces)
def __reduce__(self):
"""
@@ -349,9 +355,10 @@ def make_new_species(self, object, label="", reactive=True, check_existing=True,
orilabel = spec.label
label = orilabel
i = 2
- while any([label in phase.names for phase in self.edge.phase_system.phases.values()]):
- label = orilabel + "-" + str(i)
- i += 1
+ if self.edge.phase_system: # None when RMS not installed
+ while any([label in phase.names for phase in self.edge.phase_system.phases.values()]):
+ label = orilabel + "-" + str(i)
+ i += 1
spec.label = label
logging.debug("Creating new species %s", spec.label)
@@ -594,7 +601,7 @@ def make_new_pdep_reaction(self, forward):
return forward
- def enlarge(self, new_object=None, react_edge=False, unimolecular_react=None, bimolecular_react=None, trimolecular_react=None):
+ def enlarge(self, new_object=None, react_edge=False, unimolecular_react=None, bimolecular_react=None, trimolecular_react=None, requires_rms=False):
"""
Enlarge a reaction model by processing the objects in the list `new_object`.
If `new_object` is a
@@ -641,7 +648,7 @@ def enlarge(self, new_object=None, react_edge=False, unimolecular_react=None, bi
# Add new species
if new_species not in self.core.species:
- reactions_moved_from_edge = self.add_species_to_core(new_species)
+ reactions_moved_from_edge = self.add_species_to_core(new_species, requires_rms=requires_rms)
else:
reactions_moved_from_edge = []
@@ -649,7 +656,7 @@ def enlarge(self, new_object=None, react_edge=False, unimolecular_react=None, bi
pdep_network, new_species = new_object
new_reactions.extend(pdep_network.explore_isomer(new_species))
- self.process_new_reactions(new_reactions, new_species, pdep_network)
+ self.process_new_reactions(new_reactions, new_species, pdep_network, requires_rms=requires_rms)
else:
raise TypeError(
@@ -673,7 +680,7 @@ def enlarge(self, new_object=None, react_edge=False, unimolecular_react=None, bi
if len(products) == 1 and products[0] == species:
new_reactions = network.explore_isomer(species)
- self.process_new_reactions(new_reactions, species, network)
+ self.process_new_reactions(new_reactions, species, network, requires_rms=requires_rms)
network.update_configurations(self)
index = 0
break
@@ -698,7 +705,7 @@ def enlarge(self, new_object=None, react_edge=False, unimolecular_react=None, bi
# Identify a core species which was used to generate the reaction
# This is only used to determine the reaction direction for processing
spc = spcTuple[0]
- self.process_new_reactions(rxnList, spc)
+ self.process_new_reactions(rxnList, spc, requires_rms=requires_rms)
################################################################
# Begin processing the new species and reactions
@@ -710,7 +717,7 @@ def enlarge(self, new_object=None, react_edge=False, unimolecular_react=None, bi
# Do thermodynamic filtering
if not np.isinf(self.thermo_tol_keep_spc_in_edge) and self.new_species_list != []:
- self.thermo_filter_species(self.new_species_list)
+ self.thermo_filter_species(self.new_species_list, requires_rms=requires_rms)
# Update unimolecular (pressure dependent) reaction networks
if self.pressure_dependence:
@@ -807,7 +814,7 @@ def clear_surface_adjustments(self):
self.new_surface_spcs_loss = set()
self.new_surface_rxns_loss = set()
- def process_new_reactions(self, new_reactions, new_species, pdep_network=None, generate_thermo=True, generate_kinetics=True):
+ def process_new_reactions(self, new_reactions, new_species, pdep_network=None, generate_thermo=True, generate_kinetics=True, requires_rms=False):
"""
Process a list of newly-generated reactions involving the new core
species or explored isomer `new_species` in network `pdep_network`.
@@ -833,12 +840,12 @@ def process_new_reactions(self, new_reactions, new_species, pdep_network=None, g
if spec not in self.core.species:
all_species_in_core = False
if spec not in self.edge.species:
- self.add_species_to_edge(spec)
+ self.add_species_to_edge(spec, requires_rms=requires_rms)
for spec in rxn.products:
if spec not in self.core.species:
all_species_in_core = False
if spec not in self.edge.species:
- self.add_species_to_edge(spec)
+ self.add_species_to_edge(spec, requires_rms=requires_rms)
isomer_atoms = sum([len(spec.molecule[0].atoms) for spec in rxn.reactants])
@@ -866,9 +873,9 @@ def process_new_reactions(self, new_reactions, new_species, pdep_network=None, g
# The reaction is not new, so it should already be in the core or edge
continue
if all_species_in_core:
- self.add_reaction_to_core(rxn)
+ self.add_reaction_to_core(rxn, requires_rms=requires_rms)
else:
- self.add_reaction_to_edge(rxn)
+ self.add_reaction_to_edge(rxn, requires_rms=requires_rms)
else:
# Add the reaction to the appropriate unimolecular reaction network
# If pdep_network is not None then that will be the network the
@@ -1111,7 +1118,7 @@ def log_enlarge_summary(
logging.info(" The model edge has {0:d} species and {1:d} reactions".format(edge_species_count, edge_reaction_count))
logging.info("")
- def add_species_to_core(self, spec):
+ def add_species_to_core(self, spec, requires_rms=False):
"""
Add a species `spec` to the reaction model core (and remove from edge if
necessary). This function also moves any reactions in the edge that gain
@@ -1129,7 +1136,7 @@ def add_species_to_core(self, spec):
if spec in self.edge.species:
# remove forbidden species from edge
logging.info("Species {0} was Forbidden and not added to Core...Removing from Edge.".format(spec))
- self.remove_species_from_edge(self.reaction_systems, spec)
+ self.remove_species_from_edge(self.reaction_systems, spec, requires_rms=requires_rms)
# remove any empty pdep networks as a result of species removal
if self.pressure_dependence:
self.remove_empty_pdep_networks()
@@ -1141,7 +1148,8 @@ def add_species_to_core(self, spec):
rxn_list = []
if spec in self.edge.species:
- self.edge.phase_system.pass_species(spec.label, self.core.phase_system)
+ if requires_rms:
+ self.edge.phase_system.pass_species(spec.label, self.core.phase_system)
# If species was in edge, remove it
logging.debug("Removing species %s from edge.", spec)
self.edge.species.remove(spec)
@@ -1161,31 +1169,26 @@ def add_species_to_core(self, spec):
# Move any identified reactions to the core
for rxn in rxn_list:
- self.add_reaction_to_core(rxn)
+ self.add_reaction_to_core(rxn, requires_rms=requires_rms)
logging.debug("Moving reaction from edge to core: %s", rxn)
- else:
- if spec.molecule[0].contains_surface_site():
- self.core.phase_system.phases["Surface"].add_species(spec, edge_phase=self.edge.phase_system.phases["Surface"])
- self.edge.phase_system.species_dict[spec.label] = spec
- self.core.phase_system.species_dict[spec.label] = spec
- else:
- self.core.phase_system.phases["Default"].add_species(spec, edge_phase=self.edge.phase_system.phases["Default"])
- self.edge.phase_system.species_dict[spec.label] = spec
- self.core.phase_system.species_dict[spec.label] = spec
+ elif requires_rms:
+ destination_phase = "Surface" if spec.molecule[0].contains_surface_site() else "Default"
+ self.core.phase_system.phases[destination_phase].add_species(spec, edge_phase=self.edge.phase_system.phases[destination_phase])
+ self.edge.phase_system.species_dict[spec.label] = spec
+ self.core.phase_system.species_dict[spec.label] = spec
return rxn_list
- def add_species_to_edge(self, spec):
+ def add_species_to_edge(self, spec, requires_rms=False):
"""
- Add a species `spec` to the reaction model edge.
+ Add a species `spec` to the reaction model edge and optionally the RMS phase.
"""
self.edge.species.append(spec)
- if spec.molecule[0].contains_surface_site():
- self.edge.phase_system.phases["Surface"].add_species(spec)
- self.edge.phase_system.species_dict[spec.label] = spec
- else:
- self.edge.phase_system.phases["Default"].add_species(spec)
- self.edge.phase_system.species_dict[spec.label] = spec
+ if not requires_rms:
+ return
+ destination_phase = "Surface" if spec.molecule[0].contains_surface_site() else "Default"
+ self.edge.phase_system.phases[destination_phase].add_species(spec)
+ self.edge.phase_system.species_dict[spec.label] = spec
def set_thermodynamic_filtering_parameters(
self, Tmax, thermo_tol_keep_spc_in_edge, min_core_size_for_prune, maximum_edge_species, reaction_systems
@@ -1209,7 +1212,7 @@ def set_thermodynamic_filtering_parameters(
self.reaction_systems = reaction_systems
self.maximum_edge_species = maximum_edge_species
- def thermo_filter_species(self, spcs):
+ def thermo_filter_species(self, spcs, requires_rms=False):
"""
checks Gibbs energy of the species in species against the
maximum allowed Gibbs energy
@@ -1224,13 +1227,13 @@ def thermo_filter_species(self, spcs):
"greater than the thermo_tol_keep_spc_in_edge of "
"{3} ".format(spc, G, Gn, self.thermo_tol_keep_spc_in_edge)
)
- self.remove_species_from_edge(self.reaction_systems, spc)
+ self.remove_species_from_edge(self.reaction_systems, spc, requires_rms=requires_rms)
# Delete any networks that became empty as a result of pruning
if self.pressure_dependence:
self.remove_empty_pdep_networks()
- def thermo_filter_down(self, maximum_edge_species, min_species_exist_iterations_for_prune=0):
+ def thermo_filter_down(self, maximum_edge_species, min_species_exist_iterations_for_prune=0, requires_rms=False):
"""
removes species from the edge based on their Gibbs energy until maximum_edge_species
is reached under the constraint that all removed species are older than
@@ -1272,7 +1275,7 @@ def thermo_filter_down(self, maximum_edge_species, min_species_exist_iterations_
logging.info(
"Removing species {0} from edge to meet maximum number of edge species, Gibbs " "number is {1}".format(spc, Gns[rInds[i]])
)
- self.remove_species_from_edge(self.reaction_systems, spc)
+ self.remove_species_from_edge(self.reaction_systems, spc, requires_rms=requires_rms)
# Delete any networks that became empty as a result of pruning
if self.pressure_dependence:
@@ -1302,7 +1305,7 @@ def remove_empty_pdep_networks(self):
del self.network_dict[source]
self.network_list.remove(network)
- def prune(self, reaction_systems, tol_keep_in_edge, tol_move_to_core, maximum_edge_species, min_species_exist_iterations_for_prune):
+ def prune(self, reaction_systems, tol_keep_in_edge, tol_move_to_core, maximum_edge_species, min_species_exist_iterations_for_prune, requires_rms=False):
"""
Remove species from the model edge based on the simulation results from
the list of `reaction_systems`.
@@ -1382,7 +1385,7 @@ def prune(self, reaction_systems, tol_keep_in_edge, tol_move_to_core, maximum_ed
for index, spec in species_to_prune[0:prune_due_to_rate_counter]:
logging.info("Pruning species %s", spec)
logging.debug(" %-56s %10.4e", spec, max_edge_species_rate_ratios[index])
- self.remove_species_from_edge(reaction_systems, spec)
+ self.remove_species_from_edge(reaction_systems, spec, requires_rms=requires_rms)
if len(species_to_prune) - prune_due_to_rate_counter > 0:
logging.info(
"Pruning %d species to obtain an edge size of %d species", len(species_to_prune) - prune_due_to_rate_counter, maximum_edge_species
@@ -1390,7 +1393,7 @@ def prune(self, reaction_systems, tol_keep_in_edge, tol_move_to_core, maximum_ed
for index, spec in species_to_prune[prune_due_to_rate_counter:]:
logging.info("Pruning species %s", spec)
logging.debug(" %-56s %10.4e", spec, max_edge_species_rate_ratios[index])
- self.remove_species_from_edge(reaction_systems, spec)
+ self.remove_species_from_edge(reaction_systems, spec, requires_rms=requires_rms)
# Delete any networks that became empty as a result of pruning
if self.pressure_dependence:
@@ -1398,7 +1401,7 @@ def prune(self, reaction_systems, tol_keep_in_edge, tol_move_to_core, maximum_ed
logging.info("")
- def remove_species_from_edge(self, reaction_systems, spec):
+ def remove_species_from_edge(self, reaction_systems, spec, requires_rms=False):
"""
Remove species `spec` from the reaction model edge.
"""
@@ -1406,11 +1409,12 @@ def remove_species_from_edge(self, reaction_systems, spec):
# remove the species
self.edge.species.remove(spec)
self.index_species_dict.pop(spec.index)
- self.edge.phase_system.remove_species(spec)
+ if requires_rms:
+ self.edge.phase_system.remove_species(spec)
# clean up species references in reaction_systems
for reaction_system in reaction_systems:
- if not isinstance(reaction_system, Reactor):
+ if not requires_rms or not isinstance(reaction_system, RMSReactor):
try:
reaction_system.species_index.pop(spec)
except KeyError:
@@ -1482,7 +1486,7 @@ def remove_species_from_edge(self, reaction_systems, spec):
self.species_cache.remove(spec)
self.species_cache.append(None)
- def add_reaction_to_core(self, rxn):
+ def add_reaction_to_core(self, rxn, requires_rms=False):
"""
Add a reaction `rxn` to the reaction model core (and remove from edge if
necessary). This function assumes `rxn` has already been checked to
@@ -1491,7 +1495,8 @@ def add_reaction_to_core(self, rxn):
"""
if rxn not in self.core.reactions:
self.core.reactions.append(rxn)
- if rxn not in self.edge.reactions:
+
+ if requires_rms and rxn not in self.edge.reactions:
# If a reaction is not in edge but is going to add to core, it is either a seed mechanism or a newly generated reaction where all reactants and products are already in core
# If the reaction is in edge, then the corresponding rms_rxn was moved from edge phase to core phase in pass_species already.
rms_species_list = self.core.phase_system.get_rms_species_list()
@@ -1508,7 +1513,7 @@ def add_reaction_to_core(self, rxn):
if rxn in self.edge.reactions:
self.edge.reactions.remove(rxn)
- def add_reaction_to_edge(self, rxn):
+ def add_reaction_to_edge(self, rxn, requires_rms=False):
"""
Add a reaction `rxn` to the reaction model edge. This function assumes
`rxn` has already been checked to ensure it is supposed to be an edge
@@ -1517,6 +1522,8 @@ def add_reaction_to_edge(self, rxn):
edge).
"""
self.edge.reactions.append(rxn)
+ if not requires_rms:
+ return
rms_species_list = self.edge.phase_system.get_rms_species_list()
species_names = self.edge.phase_system.get_species_names()
bits = np.array([spc.molecule[0].contains_surface_site() for spc in rxn.reactants + rxn.products])
@@ -1573,7 +1580,7 @@ def get_stoichiometry_matrix(self):
stoichiometry[i, j] = nu
return stoichiometry.tocsr()
- def add_seed_mechanism_to_core(self, seed_mechanism, react=False):
+ def add_seed_mechanism_to_core(self, seed_mechanism, react=False, requires_rms=False):
"""
Add all species and reactions from `seed_mechanism`, a
:class:`KineticsPrimaryDatabase` object, to the model core. If `react`
@@ -1645,9 +1652,9 @@ def add_seed_mechanism_to_core(self, seed_mechanism, react=False):
# This unimolecular library reaction is flagged as `elementary_high_p` and has Arrhenius type kinetics.
# We should calculate a pressure-dependent rate for it
if len(rxn.reactants) == 1:
- self.process_new_reactions(new_reactions=[rxn], new_species=rxn.reactants[0])
+ self.process_new_reactions(new_reactions=[rxn], new_species=rxn.reactants[0], requires_rms=requires_rms)
else:
- self.process_new_reactions(new_reactions=[rxn], new_species=rxn.products[0])
+ self.process_new_reactions(new_reactions=[rxn], new_species=rxn.products[0], requires_rms=requires_rms)
# Perform species constraints and forbidden species checks
@@ -1684,7 +1691,7 @@ def add_seed_mechanism_to_core(self, seed_mechanism, react=False):
spec.get_liquid_volumetric_mass_transfer_coefficient_data()
spec.get_henry_law_constant_data()
- self.add_species_to_core(spec)
+ self.add_species_to_core(spec, requires_rms=requires_rms)
for rxn in self.new_reaction_list:
if self.pressure_dependence and rxn.is_unimolecular():
@@ -1695,7 +1702,7 @@ def add_seed_mechanism_to_core(self, seed_mechanism, react=False):
submit(spec, self.solvent_name)
rxn.fix_barrier_height(force_positive=True, solvent=self.solvent_name)
- self.add_reaction_to_core(rxn)
+ self.add_reaction_to_core(rxn, requires_rms=requires_rms)
# Check we didn't introduce unmarked duplicates
self.mark_chemkin_duplicates()
@@ -1707,7 +1714,7 @@ def add_seed_mechanism_to_core(self, seed_mechanism, react=False):
new_edge_reactions=[],
)
- def add_reaction_library_to_edge(self, reaction_library):
+ def add_reaction_library_to_edge(self, reaction_library, requires_rms=False):
"""
Add all species and reactions from `reaction_library`, a
:class:`KineticsPrimaryDatabase` object, to the model edge.
@@ -1772,9 +1779,9 @@ def add_reaction_library_to_edge(self, reaction_library):
# This unimolecular library reaction is flagged as `elementary_high_p` and has Arrhenius type kinetics.
# We should calculate a pressure-dependent rate for it
if len(rxn.reactants) == 1:
- self.process_new_reactions(new_reactions=[rxn], new_species=rxn.reactants[0])
+ self.process_new_reactions(new_reactions=[rxn], new_species=rxn.reactants[0], requires_rms=requires_rms)
else:
- self.process_new_reactions(new_reactions=[rxn], new_species=rxn.products[0])
+ self.process_new_reactions(new_reactions=[rxn], new_species=rxn.products[0], requires_rms=requires_rms)
# Perform species constraints and forbidden species checks
for spec in self.new_species_list:
@@ -1811,7 +1818,7 @@ def add_reaction_library_to_edge(self, reaction_library):
spec.get_liquid_volumetric_mass_transfer_coefficient_data()
spec.get_henry_law_constant_data()
- self.add_species_to_edge(spec)
+ self.add_species_to_edge(spec, requires_rms=requires_rms)
for rxn in self.new_reaction_list:
if not (
@@ -1826,7 +1833,7 @@ def add_reaction_library_to_edge(self, reaction_library):
)
):
# Don't add to the edge library reactions that were already processed
- self.add_reaction_to_edge(rxn)
+ self.add_reaction_to_edge(rxn, requires_rms=requires_rms)
if self.save_edge_species:
from rmgpy.chemkin import mark_duplicate_reaction
diff --git a/rmgpy/rmg/pdep.py b/rmgpy/rmg/pdep.py
index 0e1398e51a..834f0f7243 100644
--- a/rmgpy/rmg/pdep.py
+++ b/rmgpy/rmg/pdep.py
@@ -736,7 +736,7 @@ def remove_products_from_reactants(self):
self.reactants.remove(prod)
self.products = self.products_cache
- def update(self, reaction_model, pdep_settings):
+ def update(self, reaction_model, pdep_settings, requires_rms=False):
"""
Regenerate the :math:`k(T,P)` values for this partial network if the
network is marked as invalid.
@@ -917,7 +917,7 @@ def update(self, reaction_model, pdep_settings):
f'from the {rxn.library} library, and was not added to the model')
break
else:
- reaction_model.add_reaction_to_core(net_reaction)
+ reaction_model.add_reaction_to_core(net_reaction, requires_rms=requires_rms)
else:
# Check whether netReaction already exists in the edge as a LibraryReaction
for rxn in reaction_model.edge.reactions:
@@ -929,7 +929,7 @@ def update(self, reaction_model, pdep_settings):
f'from the {rxn.library} library, and was not added to the model')
break
else:
- reaction_model.add_reaction_to_edge(net_reaction)
+ reaction_model.add_reaction_to_edge(net_reaction, requires_rms=requires_rms)
# Set/update the net reaction kinetics using interpolation model
kdata = K[:, :, i, j].copy()
diff --git a/rmgpy/rmg/reactors.py b/rmgpy/rmg/reactionmechanismsimulator_reactors.py
similarity index 78%
rename from rmgpy/rmg/reactors.py
rename to rmgpy/rmg/reactionmechanismsimulator_reactors.py
index 9386fa2a20..f40ba703d6 100644
--- a/rmgpy/rmg/reactors.py
+++ b/rmgpy/rmg/reactionmechanismsimulator_reactors.py
@@ -30,48 +30,94 @@
"""
Contains classes for building RMS simulations
"""
-import numpy as np
-import sys
-import logging
import itertools
-import rmgpy.constants as constants
+import logging
+import sys
-if __debug__:
- try:
- from os.path import dirname, abspath, join, exists
+import numpy as np
- path_rms = dirname(dirname(dirname(abspath(__file__))))
- from julia.api import Julia
+import rmgpy.constants as constants
- jl = Julia(sysimage=join(path_rms, "rms.so")) if exists(join(path_rms, "rms.so")) else Julia(compiled_modules=False)
- from pyrms import rms
- from diffeqpy import de
- from julia import Main
- except Exception as e:
- import warnings
-
- warnings.warn("Unable to import Julia dependencies, original error: " + str(e), RuntimeWarning)
-else:
- from pyrms import rms
- from diffeqpy import de
- from julia import Main
+NO_JULIA = True
+try:
+ import juliacall
+ from juliacall import Main
+ Main.seval("using PythonCall")
+ Main.seval("using ReactionMechanismSimulator")
+ Main.seval("using ReactionMechanismSimulator.Sundials")
+ NO_JULIA = False
+except:
+ logging.warning("Julia import failed, RMS reactors not available.")
from rmgpy import constants
-from rmgpy.species import Species
+from rmgpy.data.kinetics.depository import DepositoryReaction
+from rmgpy.data.kinetics.family import TemplateReaction
+from rmgpy.data.solvation import SolventData
+from rmgpy.kinetics.arrhenius import (
+ Arrhenius,
+ ArrheniusBM,
+ ArrheniusChargeTransfer,
+ ArrheniusEP,
+ Marcus,
+ MultiArrhenius,
+ MultiPDepArrhenius,
+ PDepArrhenius,
+)
+from rmgpy.kinetics.chebyshev import Chebyshev
+from rmgpy.kinetics.falloff import Lindemann, ThirdBody, Troe
+from rmgpy.kinetics.kineticsdata import KineticsData
+from rmgpy.kinetics.surface import StickingCoefficient, SurfaceChargeTransfer
from rmgpy.molecule.fragment import Fragment
from rmgpy.reaction import Reaction
-from rmgpy.thermo.nasa import NASAPolynomial, NASA
-from rmgpy.thermo.wilhoit import Wilhoit
+from rmgpy.solver.termination import (
+ TerminationConversion,
+ TerminationRateRatio,
+ TerminationTime,
+)
+from rmgpy.species import Species
+from rmgpy.thermo.nasa import NASA, NASAPolynomial
from rmgpy.thermo.thermodata import ThermoData
-from rmgpy.kinetics.arrhenius import Arrhenius, ArrheniusEP, ArrheniusBM, PDepArrhenius, MultiArrhenius, MultiPDepArrhenius, ArrheniusChargeTransfer, Marcus
-from rmgpy.kinetics.kineticsdata import KineticsData
-from rmgpy.kinetics.falloff import Troe, ThirdBody, Lindemann
-from rmgpy.kinetics.chebyshev import Chebyshev
-from rmgpy.data.solvation import SolventData
-from rmgpy.kinetics.surface import StickingCoefficient, SurfaceChargeTransfer
-from rmgpy.solver.termination import TerminationTime, TerminationConversion, TerminationRateRatio
-from rmgpy.data.kinetics.family import TemplateReaction
-from rmgpy.data.kinetics.depository import DepositoryReaction
+from rmgpy.thermo.wilhoit import Wilhoit
+
+
+def to_julia(obj):
+ """
+ Convert native Python object to julia object. If the object is a Python dict, it will be converted to a Julia Dict. If the object is a Python list, it will be converted to a Julia Vector. If the object is a 1-d numpy array, it will be converted to a Julia Array. If the object is a n-d (n > 1) numpy array, it will be converted to a Julia Matrix.
+ Otherwise, the object will be returned as is. Doesn't support RMG objects.
+
+ Parameters
+ ----------
+ obj : dict | list | np.ndarray | object
+ The native Python object to convert
+
+ Returns
+ -------
+ object : Main.Dict | Main.Vector | Main.Matrix | object
+ The Julia object
+ """
+ if isinstance(obj, dict):
+ return Main.PythonCall.pyconvert(Main.Dict, obj)
+ elif isinstance(obj, (list, np.ndarray)):
+ if getattr(obj, "shape", False) and len(obj.shape) > 1:
+ return Main.PythonCall.pyconvert(Main.Matrix, obj)
+ return Main.PythonCall.pyconvert(Main.Vector, obj)
+ else: # Other native Python project does not need special conversion.
+ return obj
+
+NO_JULIA = False
+try:
+ if __debug__:
+ from os.path import abspath, dirname, exists, join
+
+ from julia.api import Julia
+ path_rms = dirname(dirname(dirname(abspath(__file__))))
+ jl = Julia(sysimage=join(path_rms, "rms.so")) if exists(join(path_rms, "rms.so")) else Julia(compiled_modules=False)
+ from diffeqpy import de
+ from julia import Main
+ from pyrms import rms
+except Exception as e:
+ logging.info("Unable to import Julia dependencies, original error: " + str(e) + ". RMS features will not be available on this execution.")
+ NO_JULIA = True
class PhaseSystem:
@@ -153,8 +199,11 @@ def pass_species(self, label, phasesys):
rxnlist = []
for i, rxn in enumerate(self.phases[phase_label].reactions):
- if (spc.name in [spec.name for spec in rxn.reactants + rxn.products]) and all(
- [spec.name in phasesys.species_dict for spec in rxn.reactants + rxn.products]
+ reactants = Main.pylist(rxn.reactants)
+ products = Main.pylist(rxn.products)
+ reacs_and_prods = reactants + products
+ if (spc.name in [spec.name for spec in reacs_and_prods]) and all(
+ [spec.name in phasesys.species_dict for spec in reacs_and_prods]
):
rxnlist.append(rxn)
@@ -168,8 +217,8 @@ def pass_species(self, label, phasesys):
for i, rxn in enumerate(interface.reactions):
if (
(spc in rxn.reactants or spc in rxn.products)
- and all([spec in phasesys.interfaces[key].species for spec in rxn.reactants])
- and all([spec in phasesys.interfaces[key].species for spec in rxn.products])
+ and all(spec.name in phasesys.species_dict for spec in rxn.reactants)
+ and all(spec.name in phasesys.species_dict for spec in rxn.products)
):
rxnlist.append(rxn)
@@ -422,7 +471,7 @@ def simulate(self, model_settings, simulator_settings, conditions):
max_edge_species_rate_ratios,
t,
x,
- ) = rms.selectobjects(
+ ) = Main.selectobjects(
core_react,
edge_react,
edge_domains,
@@ -448,7 +497,7 @@ def simulate(self, model_settings, simulator_settings, conditions):
model_settings.tol_rxn_to_core_deadend_radical,
atol=simulator_settings.atol,
rtol=simulator_settings.rtol,
- solver=de.CVODE_BDF(),
+ solver=Main.Sundials.CVODE_BDF(),
)
return (
@@ -473,9 +522,9 @@ def generate_reactor(self, phase_system):
Setup an RMS simulation for EdgeAnalysis
"""
phase = phase_system.phases["Default"]
- ig = rms.IdealGas(phase.species, phase.reactions)
- domain, y0, p = rms.ConstantVDomain(phase=ig, initialconds=self.initial_conditions)
- react = rms.Reactor(domain, y0, (0.0, self.tf), p=p)
+ ig = Main.IdealGas(phase.species, phase.reactions)
+ domain, y0, p = Main.ConstantVDomain(phase=ig, initialconds=to_julia(self.initial_conditions))
+ react = Main.Reactor(domain, y0, (0.0, self.tf), p=p)
return react, domain, [], p
@@ -491,31 +540,33 @@ def generate_reactor(self, phase_system):
surf = phase_system.phases["Surface"]
interface = list(phase_system.interfaces.values())[0]
if "mu" in self.initial_conditions["liquid"].keys():
- solv = rms.Solvent("solvent",rms.ConstantViscosity(self.initial_conditions["liquid"]["mu"]))
+ solv = Main.Solvent("solvent", Main.ConstantViscosity(self.initial_conditions["liquid"]["mu"]))
liq_initial_cond = self.initial_conditions["liquid"].copy()
del liq_initial_cond["mu"]
else:
solv = liq.solvent
liq_initial_cond = self.initial_conditions["liquid"]
- liq = rms.IdealDiluteSolution(liq.species, liq.reactions, solv, name="liquid",diffusionlimited=True)
- surf = rms.IdealSurface(surf.species, surf.reactions, surf.site_density, name="surface")
+ liq = Main.IdealDiluteSolution(liq.species, liq.reactions, solv, name="liquid",diffusionlimited=True)
+ surf = Main.IdealSurface(surf.species, surf.reactions, surf.site_density, name="surface")
liq_constant_species = [cspc for cspc in self.const_spc_names if cspc in [spc.name for spc in liq.species]]
cat_constant_species = [cspc for cspc in self.const_spc_names if cspc in [spc.name for spc in surf.species]]
- domainliq,y0liq,pliq = rms.ConstantTVDomain(phase=liq,initialconds=liq_initial_cond,constantspecies=liq_constant_species)
- domaincat,y0cat,pcat = rms.ConstantTAPhiDomain(phase=surf,initialconds=self.initial_conditions["surface"],constantspecies=cat_constant_species)
+ domainliq, y0liq, pliq = Main.ConstantTVDomain(phase=liq, initialconds=to_julia(liq_initial_cond), constantspecies=to_julia(liq_constant_species))
+ domaincat, y0cat, pcat = Main.ConstantTAPhiDomain(
+ phase=surf, initialconds=to_julia(liq_initial_cond), constantspecies=to_julia(cat_constant_species),
+ )
if interface.reactions == []:
- inter, pinter = rms.ReactiveInternalInterfaceConstantTPhi(
+ inter, pinter = Main.ReactiveInternalInterfaceConstantTPhi(
domainliq,
domaincat,
- Main.eval("using ReactionMechanismSimulator; Vector{ElementaryReaction}()"),
+ Main.seval("using ReactionMechanismSimulator; Vector{ElementaryReaction}()"),
self.initial_conditions["liquid"]["T"],
self.initial_conditions["surface"]["A"],
)
else:
- inter, pinter = rms.ReactiveInternalInterfaceConstantTPhi(
+ inter, pinter = Main.ReactiveInternalInterfaceConstantTPhi(
domainliq, domaincat, interface.reactions, self.initial_conditions["liquid"]["T"], self.initial_conditions["surface"]["A"]
)
- react, y0, p = rms.Reactor((domainliq, domaincat), (y0liq, y0cat), (0.0, self.tf), (inter,), (pliq, pcat, pinter))
+ react, y0, p = Main.Reactor((domainliq, domaincat), (y0liq, y0cat), (0.0, self.tf), (inter,), (pliq, pcat, pinter))
return react, (domainliq, domaincat), (inter,), p
@@ -542,27 +593,27 @@ def generate_reactor(self, phase_system):
Setup an RMS simulation for EdgeAnalysis
"""
phase = phase_system.phases["Default"]
- liq = rms.IdealDiluteSolution(phase.species, phase.reactions, phase.solvent)
- domain, y0, p = rms.ConstantTVDomain(phase=liq, initialconds=self.initial_conditions, constantspecies=self.const_spc_names)
+ liq = Main.IdealDiluteSolution(phase.species, phase.reactions, phase.solvent)
+ domain, y0, p = Main.ConstantTVDomain(phase=liq, initialconds=to_julia(self.initial_conditions), constantspecies=to_julia(self.const_spc_names))
interfaces = []
if self.inlet_conditions:
inlet_conditions = {key: value for (key, value) in self.inlet_conditions.items() if key != "F"}
total_molar_flow_rate = self.inlet_conditions["F"]
- inlet = rms.Inlet(domain, inlet_conditions, Main.eval("x->" + str(total_molar_flow_rate)))
+ inlet = Main.Inlet(domain, to_julia(inlet_conditions), Main.seval("x->" + str(total_molar_flow_rate)))
interfaces.append(inlet)
if self.outlet_conditions:
total_volumetric_flow_rate = self.outlet_conditions["Vout"]
- outlet = rms.VolumetricFlowRateOutlet(domain, Main.eval("x->" + str(total_volumetric_flow_rate)))
+ outlet = Main.VolumetricFlowRateOutlet(domain, Main.seval("x->" + str(total_volumetric_flow_rate)))
interfaces.append(outlet)
if self.evap_cond_conditions:
- kLA_kH_evap_cond = rms.kLAkHCondensationEvaporationWithReservoir(domain, self.evap_cond_conditions)
+ kLA_kH_evap_cond = Main.kLAkHCondensationEvaporationWithReservoir(domain, self.evap_cond_conditions)
interfaces.append(kLA_kH_evap_cond)
- react = rms.Reactor(domain, y0, (0.0, self.tf), interfaces, p=p)
+ react = Main.Reactor(domain, y0, (0.0, self.tf), interfaces, p=p)
return react, domain, interfaces, p
@@ -577,9 +628,9 @@ def generate_reactor(self, phase_system):
Setup an RMS simulation for EdgeAnalysis
"""
phase = phase_system.phases["Default"]
- ig = rms.IdealGas(phase.species, phase.reactions)
- domain, y0, p = rms.ConstantTPDomain(phase=ig, initialconds=self.initial_conditions)
- react = rms.Reactor(domain, y0, (0.0, self.tf), p=p)
+ ig = Main.IdealGas(phase.species, phase.reactions)
+ domain, y0, p = Main.ConstantTPDomain(phase=ig, initialconds=to_julia(self.initial_conditions))
+ react = Main.Reactor(domain, y0, (0.0, self.tf), p=p)
return react, domain, [], p
@@ -596,7 +647,7 @@ def to_rms(obj, species_names=None, rms_species_list=None, rmg_species=None):
A = obj._A.value_si
n = obj._n.value_si
Ea = obj._Ea.value_si
- return rms.Arrhenius(A, n, Ea, rms.EmptyRateUncertainty())
+ return Main.Arrhenius(A, n, Ea, Main.EmptyRateUncertainty())
elif isinstance(obj, ArrheniusChargeTransfer):
A = obj._A.value_si
if obj._T0.value_si != 1.0:
@@ -605,7 +656,7 @@ def to_rms(obj, species_names=None, rms_species_list=None, rmg_species=None):
Ea = obj._Ea.value_si
q = obj._alpha.value_si*obj._electrons.value_si
V0 = obj._V0.value_si
- return rms.Arrheniusq(A, n, Ea, q, V0, rms.EmptyRateUncertainty())
+ return Main.Arrheniusq(A, n, Ea, q, V0, Main.EmptyRateUncertainty())
elif isinstance(obj, SurfaceChargeTransfer):
A = obj._A.value_si
if obj._T0.value_si != 1.0:
@@ -614,28 +665,28 @@ def to_rms(obj, species_names=None, rms_species_list=None, rmg_species=None):
Ea = obj._Ea.value_si
q = obj._alpha.value_si*obj._electrons.value_si
V0 = obj._V0.value_si
- return rms.Arrheniusq(A, n, Ea, q, V0, rms.EmptyRateUncertainty())
+ return Main.Arrheniusq(A, n, Ea, q, V0, Main.EmptyRateUncertainty())
elif isinstance(obj, Marcus):
A = obj._A.value_si
n = obj._n.value_si
- return rms.Marcus(A,n,obj._lmbd_i_coefs.value_si,obj._lmbd_o.value_si, obj._wr.value_si, obj._wp.value_si, obj._beta.value_si, rms.EmptyRateUncertainty())
+ return Main.Marcus(A,n,obj._lmbd_i_coefs.value_si,obj._lmbd_o.value_si, obj._wr.value_si, obj._wp.value_si, obj._beta.value_si, Main.EmptyRateUncertainty())
elif isinstance(obj, PDepArrhenius):
- Ps = obj._pressures.value_si
- arrs = [to_rms(arr) for arr in obj.arrhenius]
- return rms.PdepArrhenius(Ps, arrs, rms.EmptyRateUncertainty())
+ Ps = to_julia(obj._pressures.value_si)
+ arrs = to_julia([to_rms(arr) for arr in obj.arrhenius])
+ return Main.PdepArrhenius(Ps, arrs, Main.EmptyRateUncertainty())
elif isinstance(obj, MultiArrhenius):
- arrs = [to_rms(arr) for arr in obj.arrhenius]
- return rms.MultiArrhenius(arrs, rms.EmptyRateUncertainty())
+ arrs = to_julia([to_rms(arr) for arr in obj.arrhenius])
+ return Main.MultiArrhenius(arrs, Main.EmptyRateUncertainty())
elif isinstance(obj, MultiPDepArrhenius):
- parrs = [to_rms(parr) for parr in obj.arrhenius]
- return rms.MultiPdepArrhenius(parrs, rms.EmptyRateUncertainty())
+ parrs = to_julia([to_rms(parr) for parr in obj.arrhenius])
+ return Main.MultiPdepArrhenius(parrs, Main.EmptyRateUncertainty())
elif isinstance(obj, Chebyshev):
Tmin = obj.Tmin.value_si
Tmax = obj.Tmax.value_si
Pmin = obj.Pmin.value_si
Pmax = obj.Pmax.value_si
- coeffs = obj.coeffs.value_si.tolist()
- return rms.Chebyshev(coeffs, Tmin, Tmax, Pmin, Pmax)
+ coeffs = to_julia(obj.coeffs.value_si)
+ return Main.Chebyshev(coeffs, Tmin, Tmax, Pmin, Pmax)
elif isinstance(obj, ThirdBody):
arrstr = arrhenius_to_julia_string(obj.arrheniusLow)
efficiencies = {rmg_species[i].label: float(val) for i, val in enumerate(obj.get_effective_collider_efficiencies(rmg_species)) if val != 1}
@@ -643,7 +694,7 @@ def to_rms(obj, species_names=None, rms_species_list=None, rmg_species=None):
for key, value in efficiencies.items():
dstr += '"' + key + '"' "=>" + str(value) + ","
dstr += "])"
- return Main.eval(
+ return Main.seval(
"using ReactionMechanismSimulator; ThirdBody(" + arrstr + ", Dict{Int64,Float64}([]), " + dstr + "," + "EmptyRateUncertainty())"
)
elif isinstance(obj, Lindemann):
@@ -654,7 +705,7 @@ def to_rms(obj, species_names=None, rms_species_list=None, rmg_species=None):
for key, value in efficiencies.items():
dstr += '"' + key + '"' "=>" + str(value) + ","
dstr += "])"
- return Main.eval(
+ return Main.seval(
"using ReactionMechanismSimulator; Lindemann("
+ arrhigh
+ ","
@@ -677,7 +728,7 @@ def to_rms(obj, species_names=None, rms_species_list=None, rmg_species=None):
for key, value in efficiencies.items():
dstr += '"' + key + '"' "=>" + str(value) + ","
dstr += "])"
- return Main.eval(
+ return Main.seval(
"using ReactionMechanismSimulator; Troe("
+ arrhigh
+ ","
@@ -703,11 +754,11 @@ def to_rms(obj, species_names=None, rms_species_list=None, rmg_species=None):
A = obj._A.value_si
n = obj._n.value_si
Ea = obj._Ea.value_si
- return rms.StickingCoefficient(A, n, Ea, rms.EmptyRateUncertainty())
+ return Main.StickingCoefficient(A, n, Ea, Main.EmptyRateUncertainty())
elif isinstance(obj, NASAPolynomial):
- return rms.NASApolynomial(obj.coeffs, obj.Tmin.value_si, obj.Tmax.value_si)
+ return Main.NASApolynomial(obj.coeffs, obj.Tmin.value_si, obj.Tmax.value_si)
elif isinstance(obj, NASA):
- return rms.NASA([to_rms(poly) for poly in obj.polynomials], rms.EmptyThermoUncertainty())
+ return Main.NASA([to_rms(poly) for poly in obj.polynomials], Main.EmptyThermoUncertainty())
elif isinstance(obj, Species):
if isinstance(obj.molecule[0], Fragment):
@@ -722,24 +773,25 @@ def to_rms(obj, species_names=None, rms_species_list=None, rmg_species=None):
atomnums[atm.element.symbol] += 1
else:
atomnums[atm.element.symbol] = 1
+ atomnums = to_julia(atomnums)
bondnum = len(mol.get_all_edges())
if not obj.molecule[0].contains_surface_site():
- rad = rms.getspeciesradius(atomnums, bondnum)
- diff = rms.StokesDiffusivity(rad)
+ rad = Main.getspeciesradius(atomnums, bondnum)
+ diff = Main.StokesDiffusivity(rad)
th = obj.get_thermo_data()
thermo = to_rms(th)
if obj.henry_law_constant_data:
- kH = rms.TemperatureDependentHenryLawConstant(Ts=obj.henry_law_constant_data.Ts, kHs=obj.henry_law_constant_data.kHs)
+ kH = Main.TemperatureDependentHenryLawConstant(Ts=to_julia(obj.henry_law_constant_data.Ts), kHs=to_julia(obj.henry_law_constant_data.kHs))
else:
- kH = rms.EmptyHenryLawConstant()
+ kH = Main.EmptyHenryLawConstant()
if obj.liquid_volumetric_mass_transfer_coefficient_data:
- kLA = rms.TemperatureDependentLiquidVolumetricMassTransferCoefficient(
- Ts=obj.liquid_volumetric_mass_transfer_coefficient_data.Ts, kLAs=obj.liquid_volumetric_mass_transfer_coefficient_data.kLAs
+ kLA = Main.TemperatureDependentLiquidVolumetricMassTransferCoefficient(
+ Ts=to_julia(obj.liquid_volumetric_mass_transfer_coefficient_data.Ts), kLAs=to_julia(obj.liquid_volumetric_mass_transfer_coefficient_data.kLAs)
)
else:
- kLA = rms.EmptyLiquidVolumetricMassTransferCoefficient()
- return rms.Species(
+ kLA = Main.EmptyLiquidVolumetricMassTransferCoefficient()
+ return Main.Species(
name=obj.label,
index=obj.index,
inchi="",
@@ -759,7 +811,7 @@ def to_rms(obj, species_names=None, rms_species_list=None, rmg_species=None):
else:
th = obj.get_thermo_data()
thermo = to_rms(th)
- return rms.Species(
+ return Main.Species(
name=obj.label,
index=obj.index,
inchi="",
@@ -768,33 +820,33 @@ def to_rms(obj, species_names=None, rms_species_list=None, rmg_species=None):
thermo=thermo,
atomnums=atomnums,
bondnum=bondnum,
- diffusion=rms.EmptyDiffusivity(),
+ diffusion=Main.EmptyDiffusivity(),
radius=0.0,
radicalelectrons=obj.molecule[0].multiplicity - 1,
molecularweight=0.0,
- henrylawconstant=rms.EmptyHenryLawConstant(),
- liquidvolumetricmasstransfercoefficient=rms.EmptyLiquidVolumetricMassTransferCoefficient(),
+ henrylawconstant=Main.EmptyHenryLawConstant(),
+ liquidvolumetricmasstransfercoefficient=Main.EmptyLiquidVolumetricMassTransferCoefficient(),
comment=obj.thermo.comment,
)
elif isinstance(obj, Reaction):
- reactantinds = [species_names.index(spc.label) for spc in obj.reactants]
- productinds = [species_names.index(spc.label) for spc in obj.products]
- reactants = [rms_species_list[i] for i in reactantinds]
- products = [rms_species_list[i] for i in productinds]
+ reactantinds = to_julia([species_names.index(spc.label) for spc in obj.reactants])
+ productinds = to_julia([species_names.index(spc.label) for spc in obj.products])
+ reactants = to_julia([rms_species_list[i] for i in reactantinds])
+ products = to_julia([rms_species_list[i] for i in productinds])
if isinstance(obj.kinetics, SurfaceChargeTransfer):
obj.set_reference_potential(300)
kinetics = to_rms(obj.kinetics, species_names=species_names, rms_species_list=rms_species_list, rmg_species=rmg_species)
radchange = sum([spc.molecule[0].multiplicity-1 for spc in obj.products]) - sum([spc.molecule[0].multiplicity-1 for spc in obj.reactants])
electronchange = -sum([spc.molecule[0].get_net_charge() for spc in obj.products]) + sum([spc.molecule[0].get_net_charge() for spc in obj.reactants])
- return rms.ElementaryReaction(index=obj.index, reactants=reactants, reactantinds=reactantinds, products=products, productinds=productinds, kinetics=kinetics, electronchange=electronchange, radicalchange=radchange, reversible=obj.reversible, pairs=[], comment=obj.kinetics.comment)
+ return Main.ElementaryReaction(index=obj.index, reactants=reactants, reactantinds=reactantinds, products=products, productinds=productinds, kinetics=kinetics, electronchange=electronchange, radicalchange=radchange, reversible=obj.reversible, pairs=to_julia([]), comment=obj.kinetics.comment)
elif isinstance(obj, SolventData):
- return rms.Solvent("solvent", rms.RiedelViscosity(float(obj.A), float(obj.B), float(obj.C), float(obj.D), float(obj.E)))
+ return Main.Solvent("solvent", Main.RiedelViscosity(float(obj.A), float(obj.B), float(obj.C), float(obj.D), float(obj.E)))
elif isinstance(obj, TerminationTime):
- return rms.TerminationTime(obj.time.value_si)
+ return Main.TerminationTime(obj.time.value_si)
elif isinstance(obj, TerminationConversion):
- return rms.TerminationConversion(to_rms(obj.species), obj.conversion)
+ return Main.TerminationConversion(to_rms(obj.species), obj.conversion)
elif isinstance(obj, TerminationRateRatio):
- return rms.TerminationRateRatio(obj.ratio)
+ return Main.TerminationRateRatio(obj.ratio)
elif isinstance(obj, tuple): # Handle TerminationConversion when the obj doesn't have thermo yet
return obj
else:
diff --git a/rmgpy/tools/fluxdiagram.py b/rmgpy/tools/fluxdiagram.py
index 8a472046e0..81ecf44f89 100644
--- a/rmgpy/tools/fluxdiagram.py
+++ b/rmgpy/tools/fluxdiagram.py
@@ -42,7 +42,7 @@
from rmgpy.kinetics.diffusionLimited import diffusion_limiter
from rmgpy.rmg.settings import SimulatorSettings
-from rmgpy.solver.base import TerminationTime, TerminationConversion
+from rmgpy.solver.base import TerminationConversion, TerminationTime
from rmgpy.solver.liquid import LiquidReactor
from rmgpy.tools.loader import load_rmg_job
@@ -308,7 +308,10 @@ def generate_flux_diagram(reaction_model, times, concentrations, reaction_rates,
'-pix_fmt', 'yuv420p', # Pixel format
'flux_diagram.avi'] # Output filename
- subprocess.check_call(command, cwd=output_directory)
+ try:
+ subprocess.run(command, check=True, capture_output=True, cwd=output_directory)
+ except subprocess.CalledProcessError as err:
+ raise RuntimeError(f"{err} {err.stderr.decode('utf8')} {err.stdout.decode('utf8')}") from err
################################################################################
diff --git a/test/arkane/arkaneFunctionalInputTest.py b/test/arkane/arkaneFunctionalInputTest.py
index 89f05f6c2e..8e1bcdb1b5 100644
--- a/test/arkane/arkaneFunctionalInputTest.py
+++ b/test/arkane/arkaneFunctionalInputTest.py
@@ -327,7 +327,7 @@ def test_load_input_file(self):
"examples",
"arkane",
"networks",
- "acetyl+O2",
+ "acetyl+O2_cse",
"input.py",
)
(
diff --git a/test/arkane/encorr/isodesmicTest.py b/test/arkane/encorr/isodesmicTest.py
index b5dc58a4b8..02e5684f44 100644
--- a/test/arkane/encorr/isodesmicTest.py
+++ b/test/arkane/encorr/isodesmicTest.py
@@ -31,7 +31,6 @@
This script contains unit tests of the :mod:`arkane.isodesmic` module.
"""
-
import numpy as np
from rmgpy.molecule import Molecule
@@ -137,12 +136,11 @@ def test_initializing_constraint_map(self):
"""
Test that the constraint map is properly initialized when a SpeciesConstraints object is initialized
"""
- caffeine_consts = SpeciesConstraints(self.caffeine, [self.butane, self.benzene])
- assert set(caffeine_consts.constraint_map.keys()) == {
- "H",
- "C",
- "O",
- "N",
+ consts = SpeciesConstraints(self.caffeine, [self.butane, self.benzene])
+ caffeine_features = consts._get_all_constraints(self.caffeine)
+ caffeine_constraint_list = [feat.__repr__() for feat in caffeine_features]
+
+ assert set(caffeine_constraint_list) == {
"C=O",
"C-N",
"C-H",
@@ -154,55 +152,69 @@ def test_initializing_constraint_map(self):
}
no_rings = SpeciesConstraints(self.caffeine, [self.butane, self.benzene], conserve_ring_size=False)
- assert set(no_rings.constraint_map.keys()) == {"H", "C", "O", "N", "C=O", "C-N", "C-H", "C=C", "C=N", "C-C"}
-
- atoms_only = SpeciesConstraints(self.caffeine, [self.butane], conserve_ring_size=False, conserve_bonds=False)
- assert set(atoms_only.constraint_map.keys()) == {"H", "C", "O", "N"}
+ caffeine_features = no_rings._get_all_constraints(self.caffeine)
+ caffeine_constraint_list = [feat.__repr__() for feat in caffeine_features]
+ assert set(caffeine_constraint_list) == {"C=O", "C-N", "C-H", "C=C", "C=N", "C-C"}
def test_enumerating_constraints(self):
"""
Test that a SpeciesConstraints object can properly enumerate the constraints of a given ErrorCancelingSpecies
"""
spcs_consts = SpeciesConstraints(self.benzene, [])
- assert set(spcs_consts.constraint_map.keys()) == {"C", "H", "C=C", "C-C", "C-H", "6_ring"}
-
- # Now that we have confirmed that the correct keys are present, overwrite the constraint map to set the order
- spcs_consts.constraint_map = {
- "H": 0,
- "C": 1,
- "C=C": 2,
- "C-C": 3,
- "C-H": 4,
- "6_ring": 5,
- }
+ benzene_features = spcs_consts._get_all_constraints(self.benzene)
+ benzene_constraint_list = [feat.__repr__() for feat in benzene_features]
+ assert set(benzene_constraint_list) == {"C=C", "C-C", "C-H", "6_ring"}
+
+ target_constraints, _ = spcs_consts._enumerate_constraints([benzene_features])
+ benzene_constraints = target_constraints
assert np.array_equal(
- spcs_consts._enumerate_constraints(self.propene),
- np.array([6, 3, 1, 1, 6, 0]),
+ benzene_constraints,
+ np.array([1, 3, 6, 3]), # ['6_ring', 'C-C', 'C-H', 'C=C']
)
+
+ spcs_consts.all_reference_species = [self.propene]
+ propene_features = spcs_consts._get_all_constraints(self.propene)
+ _, reference_constraints = spcs_consts._enumerate_constraints([benzene_features, propene_features])
+ propene_constraints = reference_constraints[0]
assert np.array_equal(
- spcs_consts._enumerate_constraints(self.butane),
- np.array([10, 4, 0, 3, 10, 0]),
+ propene_constraints,
+ np.array([0, 1, 6, 1]), # ['6_ring', 'C-C', 'C-H', 'C=C']
)
+
+ spcs_consts.all_reference_species = [self.butane]
+ butane_features = spcs_consts._get_all_constraints(self.butane)
+ _, reference_constraints = spcs_consts._enumerate_constraints([benzene_features, butane_features])
+ butane_constraints = reference_constraints[0]
assert np.array_equal(
- spcs_consts._enumerate_constraints(self.benzene),
- np.array([6, 6, 3, 3, 6, 1]),
+ butane_constraints,
+ np.array([0, 3, 10, 0]), # ['6_ring', 'C-C', 'C-H', 'C=C']
)
- # Caffeine and ethyne should return None since they have features not found in benzene
- assert spcs_consts._enumerate_constraints(self.caffeine) is None
- assert spcs_consts._enumerate_constraints(self.ethyne) is None
+ # Caffeine and ethyne should return empty list since they have features not found in benzene
+ spcs_consts.all_reference_species = [self.caffeine]
+ caffeine_features = spcs_consts._get_all_constraints(self.caffeine)
+ _, reference_constraints = spcs_consts._enumerate_constraints([benzene_features, caffeine_features])
+ assert len(reference_constraints) == 0
+
+ spcs_consts.all_reference_species = [self.ethyne]
+ ethyne_features = spcs_consts._get_all_constraints(self.ethyne)
+ _, reference_constraints = spcs_consts._enumerate_constraints([benzene_features, ethyne_features])
+ assert len(reference_constraints) == 0
def test_calculating_constraints(self):
"""
Test that a SpeciesConstraints object can properly return the target constraint vector and the constraint matrix
"""
spcs_consts = SpeciesConstraints(self.caffeine, [self.propene, self.butane, self.benzene, self.ethyne])
- assert set(spcs_consts.constraint_map.keys()) == {
- "H",
- "C",
- "O",
- "N",
+ caffeine_features = spcs_consts._get_all_constraints(self.caffeine)
+ propene_features = spcs_consts._get_all_constraints(self.propene)
+ butane_features = spcs_consts._get_all_constraints(self.butane)
+ benzene_features = spcs_consts._get_all_constraints(self.benzene)
+ ethyne_features = spcs_consts._get_all_constraints(self.ethyne)
+
+ caffeine_feature_list = [feat.__repr__() for feat in caffeine_features]
+ assert set(caffeine_feature_list) == {
"C=O",
"C-N",
"C-H",
@@ -213,36 +225,20 @@ def test_calculating_constraints(self):
"6_ring",
}
- # Now that we have confirmed that the correct keys are present, overwrite the constraint map to set the order
- spcs_consts.constraint_map = {
- "H": 0,
- "C": 1,
- "O": 2,
- "N": 3,
- "C=O": 4,
- "C-N": 5,
- "C-H": 6,
- "C=C": 7,
- "C=N": 8,
- "C-C": 9,
- "5_ring": 10,
- "6_ring": 11,
- }
-
target_consts, consts_matrix = spcs_consts.calculate_constraints()
# First, test that ethyne is not included in the reference set
assert spcs_consts.reference_species == [self.propene, self.butane, self.benzene]
# Then, test the output of the calculation
- assert np.array_equal(target_consts, np.array([10, 8, 2, 4, 2, 10, 10, 1, 1, 1, 1, 1]))
+ assert np.array_equal(target_consts, np.array([ 1, 1, 1, 10, 10, 1, 1, 2, 0, 8, 10, 4, 2])) # ['5_ring', '6_ring', 'C-C', 'C-H', 'C-N', 'C=C', 'C=N', 'C=O']
assert np.array_equal(
consts_matrix,
np.array(
[
- [6, 3, 0, 0, 0, 0, 6, 1, 0, 1, 0, 0],
- [10, 4, 0, 0, 0, 0, 10, 0, 0, 3, 0, 0],
- [6, 6, 0, 0, 0, 0, 6, 3, 0, 3, 0, 1],
+ [ 0, 0, 1, 6, 0, 1, 0, 0, 0, 3, 6, 0, 0], # ['5_ring', '6_ring', 'C-C', 'C-H', 'C-N', 'C=C', 'C=N', 'C=O']
+ [ 0, 0, 3, 10, 0, 0, 0, 0, 0, 4, 10, 0, 0], # ['5_ring', '6_ring', 'C-C', 'C-H', 'C-N', 'C=C', 'C=N', 'C=O']
+ [ 0, 1, 3, 6, 0, 3, 0, 0, 0, 6, 6, 0, 0], # ['5_ring', '6_ring', 'C-C', 'C-H', 'C-N', 'C=C', 'C=N', 'C=O']
]
),
)
@@ -255,12 +251,6 @@ class TestErrorCancelingScheme:
@classmethod
def setup_class(cls):
- try:
- import pyomo as pyo
- except ImportError:
- pyo = None
- cls.pyo = pyo
-
lot = LevelOfTheory("test")
cls.propene = ErrorCancelingSpecies(Molecule(smiles="CC=C"), (100, "kJ/mol"), lot, (105, "kJ/mol"))
cls.propane = ErrorCancelingSpecies(Molecule(smiles="CCC"), (75, "kJ/mol"), lot, (80, "kJ/mol"))
@@ -276,10 +266,12 @@ def setup_class(cls):
def test_creating_error_canceling_schemes(self):
scheme = ErrorCancelingScheme(
- self.propene,
- [self.butane, self.benzene, self.caffeine, self.ethyne],
- True,
- True,
+ target=self.propene,
+ reference_set=[self.butane, self.benzene, self.caffeine, self.ethyne],
+ isodesmic_class="rc2",
+ conserve_ring_size=True,
+ limit_charges=True,
+ limit_scope=True,
)
assert scheme.reference_species == [self.butane]
@@ -298,17 +290,11 @@ def test_find_error_canceling_reaction(self):
)
# Note that caffeine and ethyne will not be allowed, so for the full set the indices are [0, 1, 2]
- rxn, _ = scheme._find_error_canceling_reaction([0, 1, 2], milp_software=["lpsolve"])
+ rxn, _ = scheme._find_error_canceling_reaction([0, 1, 2])
assert rxn.species[self.butane] == -1
assert rxn.species[self.propane] == 1
assert rxn.species[self.butene] == 1
- if self.pyo is not None:
- rxn, _ = scheme._find_error_canceling_reaction([0, 1, 2], milp_software=["pyomo"])
- assert rxn.species[self.butane] == -1
- assert rxn.species[self.propane] == 1
- assert rxn.species[self.butene] == 1
-
def test_multiple_error_canceling_reactions(self):
"""
Test that multiple error canceling reactions can be found
@@ -328,20 +314,13 @@ def test_multiple_error_canceling_reactions(self):
)
reaction_list = scheme.multiple_error_canceling_reaction_search(n_reactions_max=20)
- assert len(reaction_list) == 20
+ assert len(reaction_list) == 6
reaction_string = reaction_list.__repr__()
# Consider both permutations of the products in the reaction string
rxn_str1 = " 1*CCC + 1*C=CCC >"
rxn_str2 = " 1*C=CCC + 1*CCC >"
assert any(rxn_string in reaction_string for rxn_string in [rxn_str1, rxn_str2])
- if self.pyo is not None:
- # pyomo is slower, so don't test as many
- reaction_list = scheme.multiple_error_canceling_reaction_search(n_reactions_max=5, milp_software=["pyomo"])
- assert len(reaction_list) == 5
- reaction_string = reaction_list.__repr__()
- assert any(rxn_string in reaction_string for rxn_string in [rxn_str1, rxn_str2])
-
def test_calculate_target_enthalpy(self):
"""
Test that ErrorCancelingScheme is able to calculate thermochemistry for the target species
@@ -360,10 +339,6 @@ def test_calculate_target_enthalpy(self):
],
)
- target_thermo, rxn_list = scheme.calculate_target_enthalpy(n_reactions_max=3, milp_software=["lpsolve"])
- assert target_thermo.value_si == 115000.0
+ target_thermo, rxn_list = scheme.calculate_target_enthalpy(n_reactions_max=3)
+ assert target_thermo.value_si == 110000.0
assert isinstance(rxn_list[0], ErrorCancelingReaction)
-
- if self.pyo is not None:
- target_thermo, _ = scheme.calculate_target_enthalpy(n_reactions_max=3, milp_software=["pyomo"])
- assert target_thermo.value_si == 115000.0
diff --git a/test/rmgpy/data/solvationTest.py b/test/rmgpy/data/solvationTest.py
index 72d0926d55..5aceb1b1eb 100644
--- a/test/rmgpy/data/solvationTest.py
+++ b/test/rmgpy/data/solvationTest.py
@@ -29,6 +29,8 @@
import os
+import pytest
+
from rmgpy import settings
from rmgpy.data.solvation import (
DatabaseError,
@@ -36,14 +38,12 @@
SolvationDatabase,
SolventLibrary,
get_critical_temperature,
- get_liquid_saturation_density,
get_gas_saturation_density,
+ get_liquid_saturation_density,
)
-from rmgpy.molecule import Molecule
-from rmgpy.rmg.main import RMG
-from rmgpy.rmg.main import Species
from rmgpy.exceptions import InputError
-import pytest
+from rmgpy.molecule import Molecule
+from rmgpy.rmg.main import RMG, Species
class TestSoluteDatabase:
@@ -141,8 +141,8 @@ def test_saturation_density(self):
"""
compound_name = "Hexane"
temp = 400 # in K
- assert round(abs(get_liquid_saturation_density(compound_name, temp) - 6383.22), 2) == 0
- assert round(abs(get_gas_saturation_density(compound_name, temp) - 162.99), 2) == 0
+ assert round(abs(get_liquid_saturation_density(compound_name, temp) - 6385.15), 2) == 0
+ assert round(abs(get_gas_saturation_density(compound_name, temp) - 163.02), 2) == 0
# Unsupported compound name
with pytest.raises(DatabaseError):
get_gas_saturation_density("Hexadecane", temp)
diff --git a/test/rmgpy/kinetics/kineticsSurfaceTest.py b/test/rmgpy/kinetics/kineticsSurfaceTest.py
index 5aeaf0b746..1f67453082 100644
--- a/test/rmgpy/kinetics/kineticsSurfaceTest.py
+++ b/test/rmgpy/kinetics/kineticsSurfaceTest.py
@@ -594,7 +594,7 @@ def test_is_temperature_valid(self):
Test the SurfaceChargeTransfer.is_temperature_valid() method.
"""
T_data = np.array([200, 400, 600, 800, 1000, 1200, 1400, 1600, 1800, 4000])
- valid_data = np.array([False, True, True, True, True, True, True, True, True, False], np.bool)
+ valid_data = np.array([False, True, True, True, True, True, True, True, True, False], bool)
for T, valid in zip(T_data, valid_data):
valid0 = self.surfchargerxn_reduction.is_temperature_valid(T)
assert valid0 == valid
diff --git a/test/rmgpy/molecule/resonanceTest.py b/test/rmgpy/molecule/resonanceTest.py
index 2130307e25..8c620ff8a3 100644
--- a/test/rmgpy/molecule/resonanceTest.py
+++ b/test/rmgpy/molecule/resonanceTest.py
@@ -31,7 +31,6 @@
from rmgpy.molecule.molecule import Molecule
from rmgpy.molecule.resonance import (
_clar_optimization,
- _clar_transformation,
generate_clar_structures,
generate_kekule_structure,
generate_optimal_aromatic_resonance_structures,
@@ -1369,15 +1368,6 @@ class ClarTest:
Contains unit tests for Clar structure methods.
"""
- def test_clar_transformation(self):
- """Test that clarTransformation generates an aromatic ring."""
- mol = Molecule().from_smiles("c1ccccc1")
- sssr = mol.get_smallest_set_of_smallest_rings()
- _clar_transformation(mol, sssr[0])
- mol.update_atomtypes()
-
- assert mol.is_aromatic()
-
def test_clar_optimization(self):
"""Test to ensure pi electrons are conserved during optimization"""
mol = Molecule().from_smiles("C1=CC=C2C=CC=CC2=C1") # Naphthalene
diff --git a/test/rmgpy/rmg/mainTest.py b/test/rmgpy/rmg/mainTest.py
index cea42c859c..2a68efc73d 100644
--- a/test/rmgpy/rmg/mainTest.py
+++ b/test/rmgpy/rmg/mainTest.py
@@ -32,15 +32,12 @@
import shutil
from unittest.mock import patch
-import pytest
-
import pandas as pd
+import pytest
-from rmgpy.rmg.main import RMG, initialize_log, make_profile_graph
-from rmgpy.rmg.main import RMG_Memory
-from rmgpy import get_path
-from rmgpy import settings
+from rmgpy import get_path, settings
from rmgpy.data.rmg import RMGDatabase
+from rmgpy.rmg.main import RMG, RMG_Memory, initialize_log, make_profile_graph
from rmgpy.rmg.model import CoreEdgeReactionModel
originalPath = get_path()
@@ -195,7 +192,10 @@ def setup_class(cls):
cls.outputDir = os.path.join(cls.testDir, "output_w_filters")
cls.databaseDirectory = settings["database.directory"]
- os.mkdir(cls.outputDir)
+ try:
+ os.mkdir(cls.outputDir)
+ except FileExistsError:
+ pass # output directory already exists
initialize_log(logging.INFO, os.path.join(cls.outputDir, "RMG.log"))
cls.rmg = RMG(
diff --git a/test/rmgpy/solver/liquidTest.py b/test/rmgpy/solver/liquidTest.py
index dd03ce323b..3559a595b3 100644
--- a/test/rmgpy/solver/liquidTest.py
+++ b/test/rmgpy/solver/liquidTest.py
@@ -29,7 +29,6 @@
import os
-
import numpy as np
from rmgpy.kinetics import Arrhenius
@@ -635,5 +634,7 @@ def teardown_class(cls):
import rmgpy.data.rmg
rmgpy.data.rmg.database = None
-
- os.remove(os.path.join(cls.file_dir, "restart_from_seed.py"))
+ try:
+ os.remove(os.path.join(cls.file_dir, "restart_from_seed.py"))
+ except FileNotFoundError:
+ pass # file will not be present if any tests failed
diff --git a/test/rmgpy/statmech/ndTorsionsTest.py b/test/rmgpy/statmech/ndTorsionsTest.py
index 80186b1f3e..5df907f074 100644
--- a/test/rmgpy/statmech/ndTorsionsTest.py
+++ b/test/rmgpy/statmech/ndTorsionsTest.py
@@ -109,4 +109,4 @@ def test_hindered_rotor_nd(self):
hdnd.read_scan()
assert round(abs(hdnd.Es[0]) - 8.58538448, 4) == 0
hdnd.fit()
- assert round(abs(hdnd.calc_partition_function(300.0)) - 2.899287634962152, 5) == 0
+ assert round(abs(hdnd.calc_partition_function(300.0)) - 2.899287634962152, 4) == 0
diff --git a/utilities.py b/utilities.py
index c8e7425818..87d39e7384 100644
--- a/utilities.py
+++ b/utilities.py
@@ -48,7 +48,6 @@ def check_dependencies():
print('{0:<15}{1:<15}{2}'.format('Package', 'Version', 'Location'))
missing = {
- 'lpsolve': _check_lpsolve(),
'openbabel': _check_openbabel(),
'pydqed': _check_pydqed(),
'pyrdl': _check_pyrdl(),
@@ -80,23 +79,6 @@ def check_dependencies():
""")
-def _check_lpsolve():
- """Check for lpsolve"""
- missing = False
-
- try:
- import lpsolve55
- except ImportError:
- print('{0:<30}{1}'.format('lpsolve55',
- 'Not found. Necessary for generating Clar structures for aromatic species.'))
- missing = True
- else:
- location = lpsolve55.__file__
- print('{0:<30}{1}'.format('lpsolve55', location))
-
- return missing
-
-
def _check_openbabel():
"""Check for OpenBabel"""
missing = False
@@ -310,7 +292,7 @@ def update_headers():
start of each file, be sure to double-check the results to make sure
important lines aren't accidentally overwritten.
"""
- shebang = """#!/usr/bin/env python-jl
+ shebang = """#!/usr/bin/env python
"""