diff --git a/Makefile b/Makefile index 3323593c68..c37fda75c2 100644 --- a/Makefile +++ b/Makefile @@ -6,13 +6,10 @@ DASPK=$(shell python -c 'import pydas.daspk; print pydas.daspk.__file__') DASSL=$(shell python -c 'import pydas.dassl; print pydas.dassl.__file__') -RDKIT_VERSION=$(shell python -c 'import rdkit; print rdkit.__version__') -.PHONY : all minimal main solver cantherm clean decython documentation QM mopac_travis +.PHONY : all minimal main solver check cantherm clean decython documentation mopac_travis -all: main solver QM - -noQM: main solver +all: main solver check minimal: python setup.py build_ext minimal --build-lib . --build-temp build --pyrex-c-in-temp @@ -39,35 +36,32 @@ endif cantherm: python setup.py build_ext cantherm --build-lib . --build-temp build --pyrex-c-in-temp -QM: - @ echo "Checking if you have symmetry..." - @ echo "symmetry -h" - @ echo "Checking you have rdkit..." - @ python -c 'import rdkit; print rdkit.__file__' - @ echo "Checking rdkit version..." -ifneq ($(RDKIT_VERSION),) - @ echo "Found rdkit version $(RDKIT_VERSION)" -else - $(error RDKit version out of date, please install RDKit version 2015.03.1 or later with InChI support); -endif - @ echo "Checking rdkit has InChI support..." - @ python -c 'from rdkit import Chem; assert Chem.inchi.INCHI_AVAILABLE, "RDKit installed without InChI Support. Please install with InChI."' +check: + @ python utilities.py check-dependencies documentation: $(MAKE) -C documentation html @ echo "Start at: documentation/build/html/index.html" clean: - python setup.py clean --build-temp build - rm -rf build/ - find . -name '*.so' -exec rm -f '{}' \; - find . -name '*.pyc' -exec rm -f '{}' \; - + @ echo "Removing build directory..." + @ python setup.py clean --build-temp build + @ echo "Removing compiled files..." + @ python utilities.py clean + @ echo "Cleanup completed." + clean-solver: - rm -r build/pyrex/rmgpy/solver/ - rm -r build/build/pyrex/rmgpy/solver/ - find rmgpy/solver/ -name '*.so' -exec rm -f '{}' \; - find rmgpy/solver/ -name '*.pyc' -exec rm -f '{}' \; + @ echo "Removing solver build directories..." +ifeq ($(OS),Windows_NT) + @ -rd /s /q build\pyrex\rmgpy\solver + @ -rd /s /q build\build\pyrex\rmgpy\solver +else + @ -rm -r build/pyrex/rmgpy/solver/ + @ -rm -r build/build/pyrex/rmgpy/solver/ +endif + @ echo "Removing compiled files..." + @ python utilities.py clean-solver + @ echo "Cleanup completed." decython: # de-cythonize all but the 'minimal'. Helpful for debugging in "pure python" mode. @@ -75,50 +69,42 @@ decython: find . -name *.pyc -exec rm -f '{}' \; test-all: -ifeq ($(OS),Windows_NT) - nosetests --nocapture --nologcapture --all-modules --verbose --with-coverage --cover-inclusive --cover-package=rmgpy --cover-erase --cover-html --cover-html-dir=testing/coverage --exe rmgpy -else +ifneq ($(OS),Windows_NT) mkdir -p testing/coverage rm -rf testing/coverage/* - nosetests --nocapture --nologcapture --all-modules --verbose --with-coverage --cover-inclusive --cover-package=rmgpy --cover-erase --cover-html --cover-html-dir=testing/coverage --exe rmgpy endif + nosetests --nocapture --nologcapture --all-modules --verbose --with-coverage --cover-inclusive --cover-package=rmgpy --cover-erase --cover-html --cover-html-dir=testing/coverage --exe rmgpy test-unittests: -ifeq ($(OS),Windows_NT) - nosetests --nocapture --nologcapture --all-modules -A 'not functional' --verbose --with-coverage --cover-inclusive --cover-package=rmgpy --cover-erase --cover-html --cover-html-dir=testing/coverage --exe rmgpy -else +ifneq ($(OS),Windows_NT) mkdir -p testing/coverage rm -rf testing/coverage/* - nosetests --nocapture --nologcapture --all-modules -A 'not functional' --verbose --with-coverage --cover-inclusive --cover-package=rmgpy --cover-erase --cover-html --cover-html-dir=testing/coverage --exe rmgpy endif + nosetests --nocapture --nologcapture --all-modules -A 'not functional' --verbose --with-coverage --cover-inclusive --cover-package=rmgpy --cover-erase --cover-html --cover-html-dir=testing/coverage --exe rmgpy test test-unittests-non-auth: -ifeq ($(OS),Windows_NT) - nosetests --nocapture --nologcapture --all-modules -A 'not functional and not auth' --verbose --with-coverage --cover-inclusive --cover-package=rmgpy --cover-erase --cover-html --cover-html-dir=testing/coverage --exe rmgpy -else +ifneq ($(OS),Windows_NT) mkdir -p testing/coverage rm -rf testing/coverage/* - nosetests --nocapture --nologcapture --all-modules -A 'not functional and not auth' --verbose --with-coverage --cover-inclusive --cover-package=rmgpy --cover-erase --cover-html --cover-html-dir=testing/coverage --exe rmgpy endif + nosetests --nocapture --nologcapture --all-modules -A 'not functional and not auth' --verbose --with-coverage --cover-inclusive --cover-package=rmgpy --cover-erase --cover-html --cover-html-dir=testing/coverage --exe rmgpy test-functional: -ifeq ($(OS),Windows_NT) - nosetests --nocapture --nologcapture --all-modules -A 'functional' --verbose --exe rmgpy -else +ifneq ($(OS),Windows_NT) mkdir -p testing/coverage rm -rf testing/coverage/* - nosetests --nocapture --nologcapture --all-modules -A 'functional' --verbose --exe rmgpy endif + nosetests --nocapture --nologcapture --all-modules -A 'functional' --verbose --exe rmgpy test-database: nosetests -v -d testing/databaseTest.py -eg0: noQM +eg0: all mkdir -p testing/eg0 rm -rf testing/eg0/* cp examples/rmg/superminimal/input.py testing/eg0/input.py @ echo "Running eg0: superminimal (H2 oxidation) example" python rmg.py testing/eg0/input.py -eg1: noQM +eg1: all mkdir -p testing/eg1 rm -rf testing/eg1/* cp examples/rmg/minimal/input.py testing/eg1/input.py @@ -168,7 +154,7 @@ eg7: all @ echo "Running eg7: gri_mech_rxn_lib example" python rmg.py testing/eg7/input.py -scoop: noQM +scoop: all mkdir -p testing/scoop rm -rf testing/scoop/* cp examples/rmg/minimal/input.py testing/scoop/input.py diff --git a/documentation/source/_static/RMG-logo.pdf b/documentation/source/_static/RMG-logo.pdf deleted file mode 100644 index 04e67a07df..0000000000 Binary files a/documentation/source/_static/RMG-logo.pdf and /dev/null differ diff --git a/documentation/source/_static/custom.css b/documentation/source/_static/custom.css new file mode 100644 index 0000000000..50044b6385 --- /dev/null +++ b/documentation/source/_static/custom.css @@ -0,0 +1,81 @@ +@import url("classic.css"); + +/* ----- sidebar styles ----- */ + +div.sidebarlogo { + text-align: center; + padding: 20px; +} + +div.sphinxsidebarwrapper { + padding: 0; +} + +#sidebar .icon { + float: left; + width: 40px; + padding: 5px; +} + +#sidebar div.menuitem { + padding: 5px; +} + +#sidebar div.menutext { + text-align: left; + margin-left: 50px; + padding: 1px 5px; +} + +#sidebar div.menutitle { + color: #F04040; + font-size: 125%; + vertical-align: bottom; +} + +#sidebar div.menudesc { + color: #FFFFFF; + font-size: 90%; + vertical-align: top; +} + +#rmgsidebar a { + text-decoration: none; +} + +#rmgsidebar div.menuitem { + min-height: 50px; +} + +#rmgsidebar div.menuitem:hover { + background-color: #606060; +} + +/* ----- table styles ----- */ + +th { + background-color: transparent; +} + +table.rmg { + border-top: 2px solid #000000; + border-bottom: 2px solid #000000; + border-collapse: collapse; + margin: 10px 10px 20px 30px; +} + +table.rmg thead { + border-bottom: 1px solid #000000; +} + +table.rmg th, +table.rmg td { + padding: 5px; +} + +/* ----- code block styles ----- */ + +pre { + border-top: 1px solid #F04040; + border-bottom: 1px solid #F04040; +} diff --git a/documentation/source/_static/database-logo-small.png b/documentation/source/_static/database-logo-small.png deleted file mode 100644 index 59e4508d41..0000000000 Binary files a/documentation/source/_static/database-logo-small.png and /dev/null differ diff --git a/documentation/source/_static/documentation-icon.png b/documentation/source/_static/documentation-icon.png deleted file mode 100644 index 88076af3ab..0000000000 Binary files a/documentation/source/_static/documentation-icon.png and /dev/null differ diff --git a/documentation/source/_static/github-icon.png b/documentation/source/_static/github-icon.png new file mode 100644 index 0000000000..119a56e22f Binary files /dev/null and b/documentation/source/_static/github-icon.png differ diff --git a/documentation/source/_static/groupdraw-logo-small.png b/documentation/source/_static/groupdraw-logo-small.png deleted file mode 100644 index 44cef2a54c..0000000000 Binary files a/documentation/source/_static/groupdraw-logo-small.png and /dev/null differ diff --git a/documentation/source/_static/input-logo-small.png b/documentation/source/_static/input-logo-small.png deleted file mode 100644 index 9542994e8e..0000000000 Binary files a/documentation/source/_static/input-logo-small.png and /dev/null differ diff --git a/documentation/source/_static/kinetics-search-logo-small.png b/documentation/source/_static/kinetics-search-logo-small.png deleted file mode 100644 index 94dd1210e9..0000000000 Binary files a/documentation/source/_static/kinetics-search-logo-small.png and /dev/null differ diff --git a/documentation/source/_static/moleculedraw-logo-small.png b/documentation/source/_static/moleculedraw-logo-small.png deleted file mode 100644 index 4bb0fc055d..0000000000 Binary files a/documentation/source/_static/moleculedraw-logo-small.png and /dev/null differ diff --git a/documentation/source/_static/pdep-logo-small.png b/documentation/source/_static/pdep-logo-small.png deleted file mode 100644 index 6297914e67..0000000000 Binary files a/documentation/source/_static/pdep-logo-small.png and /dev/null differ diff --git a/documentation/source/_static/rmg-icon.png b/documentation/source/_static/rmg-icon.png new file mode 100644 index 0000000000..231a23b689 Binary files /dev/null and b/documentation/source/_static/rmg-icon.png differ diff --git a/documentation/source/_static/rmg-logo-big.pdf b/documentation/source/_static/rmg-logo-big.pdf new file mode 100755 index 0000000000..2023eedafc Binary files /dev/null and b/documentation/source/_static/rmg-logo-big.pdf differ diff --git a/documentation/source/_static/rmg-logo-small.png b/documentation/source/_static/rmg-logo-small.png index 7e6d7742c9..741de8ed83 100644 Binary files a/documentation/source/_static/rmg-logo-small.png and b/documentation/source/_static/rmg-logo-small.png differ diff --git a/documentation/source/_static/rmg.css b/documentation/source/_static/rmg.css deleted file mode 100644 index f538f4b28e..0000000000 --- a/documentation/source/_static/rmg.css +++ /dev/null @@ -1,318 +0,0 @@ -@import url("basic.css"); - -/* -- page layout ----------------------------------------------------------- */ - -body { - font-family: sans-serif; - font-size: 100%; - background-color: #F0F0F0; - color: #000; - margin: 0; - padding: 0; -} - -div.header { - color: #993333; - font-size: 0.8em; - text-align: right; - padding: 4px 10px 4px 0; -} - -div.document { - -} - -div.documentwrapper { - float: left; - width: 100%; -} - -div.bodywrapper { - background-color: #FFFFFF; - border: 1px solid #993333; - border-right: 0px; - margin: 0 0 0 230px; -} - -div.body { - color: #000000; - padding: 0px 10px 10px 10px; -} - -div.footer { - color: #808080; - width: 100%; - padding: 4px 0 4px 0; - margin: 40px 0 0 0; - text-align: center; - font-size: 0.7em; -} - -div.footer a { - color: #808080; - text-decoration: underline; -} - -div.related { - color: #993333; - font-size: 0.8em; - padding: 4px 0 2px 0; - border-bottom: 1px solid #993333; -} - -div.related a { - color: #993333; -} - -div.related ul { - padding: 0 0 0 0; -} - -div.related ul li { - margin-right: 0px; - padding: 0 0 0 0; -} - -div.sphinxsidebar { -} - -div.sidebarlogo { - text-align: center; - margin: 0 0 20px 0; -} - -div.sphinxsidebar h3 { - font-family: sans-serif; - color: #993333; - font-size: 1.4em; - font-weight: normal; - margin: 0; - padding: 0; -} - -div.sphinxsidebar h3 a { - color: #993333; -} - -div.sphinxsidebar h4 { - font-family: sans-serif; - color: #993333; - font-size: 1.3em; - font-weight: normal; - margin: 5px 0 0 0; - padding: 0; -} - -div.sphinxsidebar p { - color: #993333; -} - -div.sphinxsidebar p.topless { - margin: 5px 10px 10px 10px; - word-wrap: break-word; -} - -div.sphinxsidebar ul { - margin: 10px; - padding: 0; - color: #ffffff; -} - -div.sphinxsidebar a { - color: #993333; -} - -div.sphinxsidebar input { - border: 1px solid #993333; - font-family: sans-serif; - font-size: 1em; -} - - - -/* -- hyperlink styles ------------------------------------------------------ */ - -a { - color: #993333; - text-decoration: none; -} - -a:visited { - color: #993333; - text-decoration: none; -} - -a:hover { - text-decoration: underline; -} - - - -/* -- body styles ----------------------------------------------------------- */ - -div.body h1, -div.body h2, -div.body h3, -div.body h4, -div.body h5, -div.body h6 { - font-family: sans-serif; - font-weight: normal; - color: #993333; - margin: 20px 0px 10px 0px; - padding: 3px 0px 3px 0px; -} - -div.body h1 { margin-top: 0; font-size: 200%; } -div.body h2 { font-size: 160%; } -div.body h3 { font-size: 140%; } -div.body h4 { font-size: 120%; } -div.body h5 { font-size: 110%; } -div.body h6 { font-size: 100%; } - -a.headerlink { - color: #c60f0f; - font-size: 0.8em; - padding: 0 4px 0 4px; - text-decoration: none; -} - -a.headerlink:hover { - background-color: #c60f0f; - color: white; -} - -div.body p, div.body dd, div.body li { - text-align: justify; - line-height: 130%; -} - -div.admonition p.admonition-title + p { - display: inline; -} - -div.admonition p { - margin-bottom: 5px; -} - -div.admonition pre { - margin-bottom: 5px; -} - -div.admonition ul, div.admonition ol { - margin-bottom: 5px; -} - -div.note { - background-color: #eee; - border: 1px solid #ccc; -} - -div.seealso { - background-color: #ffc; - border: 1px solid #ff6; -} - -div.topic { - background-color: #eee; -} - -div.warning { - background-color: #ffe4e4; - border: 1px solid #f66; -} - -p.admonition-title { - display: inline; -} - -p.admonition-title:after { - content: ":"; -} - -pre { - padding: 5px; - background-color: #eeffcc; - color: #333333; - line-height: 120%; - border: 1px solid #ac9; - border-left: none; - border-right: none; -} - -tt { - background-color: #ecf0f3; - padding: 0 1px 0 1px; - font-size: 0.95em; -} - -.warning tt { - background: #efc2c2; -} - -.note tt { - background: #d6d6d6; -} - -.viewcode-back { - font-family: sans-serif; -} - -div.viewcode-block:target { - background-color: #f4debf; - border-top: 1px solid #ac9; - border-bottom: 1px solid #ac9; -} - -p.caption { - font-size: 0.9em; -} - - -/* -- table styles ----------------------------------------------------------- */ -/* -table, -table colgroup, -table th, -table tr, -table td { - border: 0px inset #000000; - padding: 0px -} -*/ - -table.data td { - padding-right: 3px -} - -table.data, table.data thead, table.data tbody { - border-top: 1px solid #000000; - border-bottom: 1px solid #000000; - padding: 20px -} - - -/* Some custom sidebar link styling */ -div.customsidebar { - border-top: 1px solid #993333; -} - -.icon_sidebar { - float: left; -} -img.icon_sidebar { max-width: none;} - -.biglink_sidebar { - font-size: 120%; - font-weight: bold; - text-align: left; - vertical-align: bottom; -} - -.linkdesc { - font-size: 90%; - text-align: left; - vertical-align: top; - font-style: italic; -} - diff --git a/documentation/source/_static/solvation-search-logo-small.png b/documentation/source/_static/solvation-search-logo-small.png deleted file mode 100644 index 6544745349..0000000000 Binary files a/documentation/source/_static/solvation-search-logo-small.png and /dev/null differ diff --git a/documentation/source/_static/tools-icon.png b/documentation/source/_static/tools-icon.png deleted file mode 100644 index 16e327dec4..0000000000 Binary files a/documentation/source/_static/tools-icon.png and /dev/null differ diff --git a/documentation/source/_templates/customsidebar.html b/documentation/source/_templates/customsidebar.html deleted file mode 100644 index b4cfc6f221..0000000000 --- a/documentation/source/_templates/customsidebar.html +++ /dev/null @@ -1,73 +0,0 @@ -

-

-

- - - - - - - - - - - - - - - - - - - - - - - - - - - -
- Documentation - - -
Learn more about the RMG software
-
- Database - - -
Browse the RMG database of chemical and kinetic parameters
-
- Create Input File - - -
Online form for making an RMG input file
-
- Molecule Search - - -
Draw a molecule from its adjlist and search or estimate its properties
-
- Kinetics Search - - -
Search or estimate kinetics of a chemical reaction
-
- Solvation Search - - -
Search or estimate solvation properties of solvent and a solute
-
- Pressure Dependent Networks - - -
CanTherm pdep kinetic calculations
-
- Simulation & Tools - - -
Additional tools to supplement RMG
-
- -

diff --git a/documentation/source/_templates/index.html b/documentation/source/_templates/index.html index 2bab1dcab9..4b220a9ecb 100644 --- a/documentation/source/_templates/index.html +++ b/documentation/source/_templates/index.html @@ -2,68 +2,71 @@ {% set title = 'Overview' %} {% block body %} -

Reaction Mechanism Generator Documentation

+

+ Reaction Mechanism Generator (RMG) is an automatic chemical reaction mechanism generator + that constructs kinetic models composed of elementary chemical reaction steps using a general understanding + of how molecules react. +

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
I want to...Resource
analyze models & search databasesrmg website resources (no download needed)
make transition state theory calculationsRun Cantherm with the Canterm User's Guide
create mechanisms automaticallyDownload RMG with the RMG User's Guide
post an issue with RMGgithub issues page
contribute to RMG projectRMG developer's wiki
- +

+ This documentation site if for the newer Python version of RMG called RMG-Py. + It is developed on GitHub as RMG-Py. +

-

- Reaction Mechanism Generator (RMG) is an automatic chemical reaction mechanism generator - that constructs kinetic models composed of elementary chemical reaction - steps using a general understanding of how molecules react. +

+ CanTherm is developed and distributed as part of RMG-Py, but can be used as a stand-alone + application for Thermochemistry, Transition State Theory, and Master Equation chemical kinetics + calculations.

- +

- These pages describe the newer Python version of RMG that we call RMG-Py. - It is developed on github as RMG-Py. + RMG is free, open-source software, and is under active and ongoing development. + If you are interested in assisting with or guiding its development, please let us know. + Any issues, concerns, or comments can be posted to our + Github issues page.

- +

- CanTherm is developed and distributed as part of RMG-Py, but can be used as a stand-alone - application for Thermochemistry, Transition State Theory, and Master Equation chemical kinetics - calculations. + While you are free to use or modify RMG components as you wish, we would like to + coordinate with you. Please contact Professor William H. Green, + Professor Richard H. West, + or the The RMG Development Team for further information.

-RMG is free, open-source software, and is under active and ongoing development. -If you are interested in assisting with or guiding its development, please let us know. Any issues, -concerns, or comments can be posted to our Github issues page. -While you are free to use/modify RMG components as you wish, we would like to -co-ordinate with you. Please contact Professor William H. Green, -Professor Richard H. West, -or the The RMG Development Team for further information. +

Quick reference guide:

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
I want to...Resource
analyze models & search databasesRMG website resources (no download needed)
make transition state theory calculationsRun Cantherm with the Canterm User's Guide
create mechanisms automaticallyDownload RMG with the RMG User's Guide
post an issue with RMGGitHub issues page
contribute to RMG projectRMG developer's wiki

Parts of the documentation:

diff --git a/documentation/source/_templates/layout.html b/documentation/source/_templates/layout.html index f9b27a9d06..eb1f46a9f7 100644 --- a/documentation/source/_templates/layout.html +++ b/documentation/source/_templates/layout.html @@ -11,39 +11,60 @@ {% endblock %} {% block rootrellink %} -
  • Documentation {{ release }} »
  • +
  • RMG »
  • +
  • {{ release }} Documentation »
  • {% endblock %} -{%- block header %} -
    -  -
    -{%- endblock %} - -{%- block sidebarlogo %} - -{%- endblock %} - -{%- block relbar1 %}{% endblock %} - - {%- block document %} -
    - {%- if render_sidebar %} -
    - {%- endif %} -
    - {{ relbar() }} - {% block body %} {% endblock %} -
    - {%- if render_sidebar %} +{%- block sidebar2 %} + {%- if render_sidebar %} + + {%- endif %} +{% endblock %} diff --git a/documentation/source/_templates/sourcelink.html b/documentation/source/_templates/sourcelink.html index f4eda6dca9..2a1ee0fe83 100644 --- a/documentation/source/_templates/sourcelink.html +++ b/documentation/source/_templates/sourcelink.html @@ -8,7 +8,12 @@ :license: BSD, see LICENSE for details. #} {%- if show_source and sourcename %} - -

    Edit this page

    +
    +

    {{ _('This Page') }}

    + +
    {%- endif %} diff --git a/documentation/source/conf.py b/documentation/source/conf.py index 37a8018bb4..5b334ba2b8 100644 --- a/documentation/source/conf.py +++ b/documentation/source/conf.py @@ -27,7 +27,7 @@ # Add any Sphinx extension module names here, as strings. They can be extensions # coming with Sphinx (named 'sphinx.ext.*') or your custom ones. -extensions = ['sphinx.ext.autodoc', 'sphinx.ext.mathjax'] +extensions = ['sphinx.ext.autodoc', 'sphinx.ext.napoleon', 'sphinx.ext.mathjax'] # Add any paths that contain templates here, relative to this directory. templates_path = ['_templates'] @@ -42,8 +42,10 @@ master_doc = 'contents' # General information about the project. +import datetime +year = datetime.datetime.now().year project = u'RMG Py' -copyright = u'2002-2017, William H. Green, Richard H. West, and the RMG Team' +copyright = u'2002-{0}, William H. Green, Richard H. West, and the RMG Team'.format(year) # The version info for the project you're documenting, acts as replacement for # |version| and |release|, also used in various other places throughout the @@ -102,20 +104,43 @@ # Theme options are theme-specific and customize the look and feel of a theme # further. For a list of options available for each theme, see the # documentation. -#html_theme_options = {} +html_theme_options = { + 'sidebarwidth': 250, + # 'collapsiblesidebar': True, + 'footerbgcolor': '#404040', # Background color for the footer line. + 'footertextcolor': '#FFFFFF', # Text color for the footer line. + 'sidebarbgcolor': '#404040', # Background color for the sidebar. + 'sidebarbtncolor': '#505050', # Background color for the sidebar collapse button (used when collapsiblesidebar is True). + 'sidebartextcolor': '#FFFFFF', # Text color for the sidebar. + 'sidebarlinkcolor': '#F04040', # Link color for the sidebar. + 'relbarbgcolor': '#303030', # Background color for the relation bar. + 'relbartextcolor': '#FFFFFF', # Text color for the relation bar. + 'relbarlinkcolor': '#F04040', # Link color for the relation bar. + 'bgcolor': '#FFFFFF', # Body background color. + 'textcolor': '#000000', # Body text color. + 'linkcolor': '#F04040', # Body link color. + 'visitedlinkcolor': '#F02020', # Body color for visited links. + 'headbgcolor': '#E0E0E0', # Background color for headings. + 'headtextcolor': '#404040', # Text color for headings. + 'headlinkcolor': '#F04040', # Link color for headings. + 'codebgcolor': '#FAF0F0', # Background color for code blocks. + 'codetextcolor': '#000000', # Default text color for code blocks, if not set differently by the highlighting style. + 'bodyfont': 'sans-serif', # Font for normal text. + 'headfont': 'sans-serif', # Font for headings. +} # stylesheet to use -html_style = 'rmg.css' +html_style = 'custom.css' # Add any paths that contain custom themes here, relative to this directory. #html_theme_path = [] # The name for this set of Sphinx documents. If None, it defaults to # " v documentation". -#html_title = None +html_title = 'RMG-Py 2.1.0 Documentation' # A shorter title for the navigation bar. Default is the same as html_title. -#html_short_title = None +html_short_title = '2.1.0 Documentation' # The name of an image file (relative to this directory) to place at the top # of the sidebar. @@ -140,7 +165,9 @@ #html_use_smartypants = True # Custom sidebar templates, maps document names to template names. -html_sidebars = {'**': ['localtoc.html', 'searchbox.html', 'sourcelink.html', 'customsidebar.html']} +html_sidebars = { + 'index': ['globaltoc.html', 'searchbox.html'], +} # Additional templates that should be rendered to pages, maps page names to # template names. @@ -179,7 +206,7 @@ # If true, documentation source files will be copied to the build directory, # and a link to view the source for the current page will be available in the # sidebar. Default is True. -html_copy_source = True +#html_copy_source = True # -- Options for LaTeX output -------------------------------------------------- @@ -200,7 +227,7 @@ # The name of an image file (relative to this directory) to place at the top of # the title page. -latex_logo = '_static/RMG-logo.pdf' +latex_logo = '_static/rmg-logo-big.pdf' # For "manual" documents, if this is true, then toplevel headings are parts, # not chapters. diff --git a/documentation/source/reference/cantherm/geometry.rst b/documentation/source/reference/cantherm/geometry.rst deleted file mode 100644 index ceadd72fab..0000000000 --- a/documentation/source/reference/cantherm/geometry.rst +++ /dev/null @@ -1,5 +0,0 @@ -******************************** -rmgpy.cantherm.geometry.Geometry -******************************** - -.. autoclass:: rmgpy.cantherm.geometry.Geometry diff --git a/documentation/source/reference/cantherm/index.rst b/documentation/source/reference/cantherm/index.rst index 86a875ee04..a611c0be46 100644 --- a/documentation/source/reference/cantherm/index.rst +++ b/documentation/source/reference/cantherm/index.rst @@ -23,18 +23,32 @@ Class Description -Geometry -======== +Reading Q-Chem log files +======================== -.. currentmodule:: rmgpy.cantherm.geometry +.. currentmodule:: rmgpy.cantherm.qchem =============================== ================================================ Class Description =============================== ================================================ -:class:`Geometry` The three-dimensional geometry of a molecular conformation +:class:`QchemLog` Extract chemical parameters from Q-Chem log files =============================== ================================================ + +Reading Molpro log files +======================== + +.. currentmodule:: rmgpy.cantherm.molpro + +=============================== ================================================ +Class Description +=============================== ================================================ +:class:`MolproLog` Extract chemical parameters from Molpro log files +=============================== ================================================ + + + Input ===== @@ -65,24 +79,12 @@ Class Description -Exceptions -========== - -.. currentmodule:: rmgpy.cantherm - -=============================== ================================================ -Exception Description -=============================== ================================================ -:exc:`GaussianError` Raised when an error occurs while working with a Gaussian log file -=============================== ================================================ - - - .. toctree:: :hidden: gaussianlog - geometry + qchemlog + molprolog input kinetics main diff --git a/documentation/source/reference/cantherm/molprolog.rst b/documentation/source/reference/cantherm/molprolog.rst new file mode 100644 index 0000000000..29c1b0bb64 --- /dev/null +++ b/documentation/source/reference/cantherm/molprolog.rst @@ -0,0 +1,5 @@ +******************************* +rmgpy.cantherm.molpro.MolproLog +******************************* + +.. autoclass:: rmgpy.cantherm.molpro.MolproLog diff --git a/documentation/source/reference/cantherm/qchemlog.rst b/documentation/source/reference/cantherm/qchemlog.rst new file mode 100644 index 0000000000..add3620bc9 --- /dev/null +++ b/documentation/source/reference/cantherm/qchemlog.rst @@ -0,0 +1,5 @@ +***************************** +rmgpy.cantherm.qchem.QchemLog +***************************** + +.. autoclass:: rmgpy.cantherm.qchem.QchemLog diff --git a/documentation/source/reference/chemkin/index.rst b/documentation/source/reference/chemkin/index.rst index d6eb22a113..05a27a9ccf 100644 --- a/documentation/source/reference/chemkin/index.rst +++ b/documentation/source/reference/chemkin/index.rst @@ -52,19 +52,6 @@ Function Description -Exceptions -========== - -.. currentmodule:: rmgpy.chemkin - -=================================== ============================================ -Exception Description -=================================== ============================================ -:exc:`ChemkinError` Raised when an error occurs while working with a Chemkin file -=================================== ============================================ - - - .. toctree:: :hidden: diff --git a/documentation/source/reference/data/index.rst b/documentation/source/reference/data/index.rst index 4636c7915c..e6070d5746 100644 --- a/documentation/source/reference/data/index.rst +++ b/documentation/source/reference/data/index.rst @@ -98,22 +98,6 @@ Class/Function Description -Exceptions -========== - -.. currentmodule:: rmgpy.data - -=================================== ============================================ -Exception Description -=================================== ============================================ -:exc:`DatabaseError` Raised when an error occurs while working with the database -:exc:`InvalidActionError` Raised when an error occurs while applying a reaction recipe -:exc:`UndeterminableKineticsError` Raised when the kinetics of a given reaction cannot be determined -:exc:`StatmechFitError` Raised when an error occurs while fitting internal degrees of freedom to heat capacity data -=================================== ============================================ - - - .. toctree:: :hidden: diff --git a/documentation/source/reference/exceptions.rst b/documentation/source/reference/exceptions.rst new file mode 100644 index 0000000000..222969036b --- /dev/null +++ b/documentation/source/reference/exceptions.rst @@ -0,0 +1,7 @@ +**************************************** +RMG Exceptions (:mod:`rmgpy.exceptions`) +**************************************** + +.. currentmodule:: rmgpy.exceptions + +.. automodule:: rmgpy.exceptions diff --git a/documentation/source/reference/index.rst b/documentation/source/reference/index.rst index e747dfa4f7..57ee3cf28c 100644 --- a/documentation/source/reference/index.rst +++ b/documentation/source/reference/index.rst @@ -26,6 +26,7 @@ Module Description :mod:`rmgpy.species` Chemical species :mod:`rmgpy.statmech` Statistical mechanics models of molecular degrees of freedom :mod:`rmgpy.thermo` Thermodynamics models of chemical species +:mod:`rmgpy.exceptions` Custom RMG exception classes ======================= ======================================================== .. toctree:: @@ -46,3 +47,4 @@ Module Description species/index statmech/index thermo/index + exceptions diff --git a/documentation/source/reference/molecule/atomtype.rst b/documentation/source/reference/molecule/atomtype.rst index e0eefc54e7..765419ad36 100644 --- a/documentation/source/reference/molecule/atomtype.rst +++ b/documentation/source/reference/molecule/atomtype.rst @@ -20,14 +20,14 @@ We define the following basic atom types: Atom type Description =============== ============================================================================================================================================================== *General atom types* ---------------- -------------------------------------------------------------------------------------------------------------------------------------------------------------- +------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ ``R`` any atom with any local bond structure ``R!H`` any non-hydrogen atom with any local bond structure *Hydrogen atom types* ---------------- -------------------------------------------------------------------------------------------------------------------------------------------------------------- +------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ ``H`` hydrogen atom with up to one single bond *Carbon atom types* ---------------- -------------------------------------------------------------------------------------------------------------------------------------------------------------- +------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ ``C`` carbon atom with any local bond structure ``Ca`` carbon atom with two lone pairs and no bonds ``Cs`` carbon atom with up to four single bonds @@ -46,7 +46,7 @@ Atom type Description ``C2dc`` charged carbon atom with one lone pair (valance 2), one double bond and up to one single bond ``C2tc`` charged carbon atom with one lone pair (valance 2), one triple bond *Nitrogen atom types* ---------------- -------------------------------------------------------------------------------------------------------------------------------------------------------------- +------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ ``N`` nitrogen atom with any local bond structure ``N0sc`` charged nitrogen atom with three lone pairs (valance 0) with up to one single bond ``N1s`` nitrogen atom with two lone pairs (valance 1) and up to one single bond @@ -65,7 +65,7 @@ Atom type Description ``N5b`` nitrogen atom with with no lone pairs (valance 5) and two benzene bonds (one of the lone pairs also participates in the aromatic bond) and up to one single bond ``N5bd`` nitrogen atom with with no lone pairs (valance 5), two benzene bonds, and one double bond *Oxygen atom types* ---------------- -------------------------------------------------------------------------------------------------------------------------------------------------------------- +------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ ``O`` oxygen atom with any local bond structure ``Oa`` oxygen atom with three lone pairs and no bonds ``O0sc`` charged oxygen with three lone pairs (valance 0) and up to one single bond @@ -78,7 +78,7 @@ Atom type Description ``O4tc`` charged oxygen atom with one one pair (valance 4) and one triple bond ``O4b`` oxygen atom with one one pair (valance 4) and and two benzene bonds *Silicon atom types* ---------------- -------------------------------------------------------------------------------------------------------------------------------------------------------------- +------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ ``Si`` silicon atom with any local bond structure ``Sis`` silicon atom with four single bonds ``Sid`` silicon atom with one double bond (to carbon) and two single bonds @@ -88,7 +88,7 @@ Atom type Description ``Sib`` silicon atom with two benzene bonds and one single bond ``Sibf`` silicon atom with three benzene bonds *Sulfur atom types* ---------------- -------------------------------------------------------------------------------------------------------------------------------------------------------------- +------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ ``S`` sulfur atom with any local bond structure ``Sa`` sulfur atom with three lone pairs and no bonds ``S0sc`` charged sulfur atom with three lone pairs (valance 0) and up to one single bonds @@ -115,4 +115,12 @@ Atom type Description ``S6td`` sulfur atom with no lone pairs (valance 6), one triple bond, one double bond and up to one single bond ``S6tt`` sulfur atom with no lone pairs (valance 6) and two triple bonds ``S6tdc`` charged sulfur atom with no lone pairs (valance 6), one to two triple bonds, up to two double bonds, and up to four single bonds +*Chlorine atom types* +------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ +``Cl`` chlorine atom with any local bond structure +``Cl1s`` chlorine atom with three lone pairs and zero to one single bonds +*Iodine atom types* +------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ +``I`` iodine atom with any local bond structure +``I1s`` iodine atom with three lone pairs and zero to one single bonds =============== ============================================================================================================================================================== diff --git a/documentation/source/reference/molecule/generator.rst b/documentation/source/reference/molecule/generator.rst new file mode 100644 index 0000000000..3f9396c8c7 --- /dev/null +++ b/documentation/source/reference/molecule/generator.rst @@ -0,0 +1,5 @@ +************************ +rmgpy.molecule.generator +************************ + +.. automodule:: rmgpy.molecule.generator diff --git a/documentation/source/reference/molecule/index.rst b/documentation/source/reference/molecule/index.rst index 9e82fa8b0e..1518f41507 100644 --- a/documentation/source/reference/molecule/index.rst +++ b/documentation/source/reference/molecule/index.rst @@ -82,6 +82,24 @@ Class Description ======================= ======================================================== + +Molecule Utilities +================== + +.. currentmodule:: rmgpy.molecule + +================================ ======================================================== +Class Description +================================ ======================================================== +:mod:`rmgpy.molecule.resonance` Resonance structure generation methods +:mod:`rmgpy.molecule.kekulize` Kekule structure generation +:mod:`rmgpy.molecule.pathfinder` Resonance path enumeration +:mod:`rmgpy.molecule.generator` Molecule string representation generator +:mod:`rmgpy.molecule.parser` Molecule string representation parser +================================ ======================================================== + + + Adjacency lists =============== @@ -126,20 +144,6 @@ Class Description ======================== ======================================================= -Exceptions -========== - -.. currentmodule:: rmgpy.molecule - -=================================== ============================================ -Exception Description -=================================== ============================================ -:exc:`ElementError` Raised when an error occurs while working with chemical elements -:exc:`AtomTypeError` Raised when an error occurs while working with atom types -:exc:`InvalidAdjacencyListError` Raised when an invalid adjacency list is encountered -:exc:`ActionError` Raised when an error occurs while working with a reaction recipe action -=================================== ============================================ - .. toctree:: :hidden: @@ -157,6 +161,11 @@ Exception Description groupatom groupbond group + resonance + kekulize + pathfinder + generator + parser adjlist symmetry moleculedrawer diff --git a/documentation/source/reference/molecule/kekulize.rst b/documentation/source/reference/molecule/kekulize.rst new file mode 100644 index 0000000000..b2aee942ca --- /dev/null +++ b/documentation/source/reference/molecule/kekulize.rst @@ -0,0 +1,5 @@ +*********************** +rmgpy.molecule.kekulize +*********************** + +.. automodule:: rmgpy.molecule.kekulize diff --git a/documentation/source/reference/molecule/parser.rst b/documentation/source/reference/molecule/parser.rst new file mode 100644 index 0000000000..4d78304e3b --- /dev/null +++ b/documentation/source/reference/molecule/parser.rst @@ -0,0 +1,5 @@ +********************* +rmgpy.molecule.parser +********************* + +.. automodule:: rmgpy.molecule.parser diff --git a/documentation/source/reference/molecule/pathfinder.rst b/documentation/source/reference/molecule/pathfinder.rst new file mode 100644 index 0000000000..75572a3a40 --- /dev/null +++ b/documentation/source/reference/molecule/pathfinder.rst @@ -0,0 +1,5 @@ +************************* +rmgpy.molecule.pathfinder +************************* + +.. automodule:: rmgpy.molecule.pathfinder diff --git a/documentation/source/reference/molecule/resonance.rst b/documentation/source/reference/molecule/resonance.rst new file mode 100644 index 0000000000..88799a3c2c --- /dev/null +++ b/documentation/source/reference/molecule/resonance.rst @@ -0,0 +1,5 @@ +************************ +rmgpy.molecule.resonance +************************ + +.. automodule:: rmgpy.molecule.resonance diff --git a/documentation/source/reference/pdep/index.rst b/documentation/source/reference/pdep/index.rst index 61564add47..e32d4d59ee 100644 --- a/documentation/source/reference/pdep/index.rst +++ b/documentation/source/reference/pdep/index.rst @@ -102,18 +102,6 @@ Function Description -Exceptions -========== - -======================================= ======================================== -Exception Description -======================================= ======================================== -:exc:`NetworkError` Raised when an error occurs while working with a pressure-dependent reaction network -:exc:`InvalidMicrocanonicalRateError` Raised when the :math:`k(E)` is not consistent with the high pressure-limit kinetics or thermodynamics -======================================= ======================================== - - - .. toctree:: :hidden: diff --git a/documentation/source/reference/reaction/index.rst b/documentation/source/reference/reaction/index.rst index b7d2a48bcf..5d4d62674d 100644 --- a/documentation/source/reference/reaction/index.rst +++ b/documentation/source/reference/reaction/index.rst @@ -22,17 +22,6 @@ Class Description -Exceptions -========== - -======================= ======================================================== -Class Description -======================= ======================================================== -:class:`ReactionError` Raised when an error occurs while working with reactions -======================= ======================================================== - - - .. toctree:: :hidden: diff --git a/documentation/source/reference/rmg/index.rst b/documentation/source/reference/rmg/index.rst index 4ca7123748..a0a6af9a58 100644 --- a/documentation/source/reference/rmg/index.rst +++ b/documentation/source/reference/rmg/index.rst @@ -78,21 +78,6 @@ Class Description -Exceptions -========== - -.. currentmodule:: rmgpy.rmg - -=============================== ================================================ -Exception Description -=============================== ================================================ -:exc:`InputError` Raised when an error occurs while working with an RMG input file -:exc:`OutputError` Raised when an error occurs while saving an RMG output file -:exc:`PressureDependenceError` Raised when an error occurs while computing pressure-dependent rate coefficients -=============================== ================================================ - - - .. toctree:: :hidden: diff --git a/documentation/source/reference/solver/index.rst b/documentation/source/reference/solver/index.rst index e43a699913..e0a7af4e0e 100644 --- a/documentation/source/reference/solver/index.rst +++ b/documentation/source/reference/solver/index.rst @@ -21,6 +21,7 @@ Class Description =========================== ==================================================== :class:`ReactionSystem` Base class for all reaction systems :class:`SimpleReactor` A simple isothermal, isobaric, well-mixed batch reactor +:class:`LiquidReactor` A homogeneous, isothermal, isobaric liquid batch reactor =========================== ==================================================== @@ -42,5 +43,6 @@ Class Description reactionsystem simplereactor + liquidreactor termination diff --git a/documentation/source/reference/solver/liquidreactor.rst b/documentation/source/reference/solver/liquidreactor.rst new file mode 100644 index 0000000000..504e302e1b --- /dev/null +++ b/documentation/source/reference/solver/liquidreactor.rst @@ -0,0 +1,5 @@ +************************** +rmgpy.solver.LiquidReactor +************************** + +.. autoclass:: rmgpy.solver.LiquidReactor diff --git a/documentation/source/reference/species/index.rst b/documentation/source/reference/species/index.rst index 9305410e07..87f6edff62 100644 --- a/documentation/source/reference/species/index.rst +++ b/documentation/source/reference/species/index.rst @@ -35,17 +35,6 @@ Class Description -Exceptions -========== - -======================= ======================================================== -Class Description -======================= ======================================================== -:class:`SpeciesError` Raised when an error occurs while working with species -======================= ======================================================== - - - .. toctree:: :hidden: diff --git a/documentation/source/theory/rmg/dynamics.rst b/documentation/source/theory/rmg/dynamics.rst new file mode 100644 index 0000000000..a2aae7fdd2 --- /dev/null +++ b/documentation/source/theory/rmg/dynamics.rst @@ -0,0 +1,97 @@ +.. _dynamics: + +******************* +Dynamics Criterion +******************* + +When dealing with more complex chemical mechanisms the standard RMG flux criterion +has trouble picking up key chain branching reactions and has limited guarantees +that it accurately represents the concentrations of all species. The dynamics +criterion is a measure of how much a given reaction affects core concentrations. +This allows it to pick up key low-flux chain branching reactions and better represent +species concentrations. + +Calculating the Dynamics Criterion +================================== +Let us define rates of production :math:`P_i(t)` and consumption :math:`L_i(t)` for a given species +:math:`\frac{dc_i}{dt} = P_i(t) - L_i(t)` + +Let us define a dimensionless concentration variable we will refer to as the +accumulation number Ac for a given species + +:math:`Ac_{spc,i} = \frac{P_i}{L_i} \approx \frac{\bar{c_{i}}}{c_{i0}}` + +where :math:`\bar{c_i}` is the steady state concentration or more specifically the +concentration at which :math:`P_i = L_i` assuming :math:`L_i` scales with :math:`c_i` +and :math:`c_{i0}` is the current concentration. + +This species accumulation number is a measure of how far species i is from steady state. + +Since this number can only be calculated for core species, by itself it is only a +measure of the behavior of species i within the reaction network. + +However if we consider models with and without some edge reaction j we can define + +:math:`\Pi_{Ac,i,j} = \frac{Ac_{spc,i,withj}}{Ac_{spc,i,withoutj}}` + +Which is a measure of how much the concentration of species i is impacted by +reaction j. + +In order to directly compare multiple reactions we can then sum over all +core species involved in reaction j to get our criterion the dynamics number. + +:math:`\sum_{i\in core} |Ln(\Pi_{Ac,i,j})| = Dy > \epsilon` + +Surface Algorithm +================= +One common issue with the dynamics criterion is that it treats all core species equally. +Because of this, if the dynamics criterion is set too low it enters a feedback loop where +it adds species and then since it can't get those species' concentrations right it adds +more species and so on. In order to avoid this feedback loop the surface algorithm was developed. +It creates a new partition called the *surface* that is considered part of the core. We will +refer to the part of the core that is not part of the surface as the *bulk core*. When +operating without the dynamics criterion everything moves from edge to the bulk core as usual; +however the dynamics criterion is managed differently. When using the surface algorithm most +reactions pulled in by the dynamics criterion enter the surface instead of the bulk core. +However, unlike movement to bulk core a constraint is placed on movement to the surface. +Any reaction moved to the surface must have either both reactants or both products +in the bulk core. This prevents the dynamics criterion from pulling in reactions +to get the concentrations of species in the surface right avoiding the feedback loop. +To avoid important species being trapped in the surface we also add criteria +for movement from surface to bulk core based on flux or dynamics criterion. +However, to avoid important species being trapped in the surface we also add criteria +for movement from surface to bulk core based on flux or dynamics criterion. + +.. figure:: surfacediagram.png + +Key Parameters for Dynamics Criterion and Surface Algorithm +=========================================================== + +* **toleranceMoveEdgeReactionToCore** + + An edge reaction will be pulled directly into the bulk core if its dynamics number + ever exceeds this value. + +* **toleranceMoveEdgeReactionToSurface** + + An edge reaction will be pulled into the surface if its dynamics number + ever exceeds this value. + +* **toleranceMoveEdgeReactionToCoreInterrupt** + + When any reaction's dynamics number exceeds this value the simulation will be interrupted. + +* **toleranceMoveEdgeReactionToSurfaceInterrupt** + + When the dynamics number of any reaction that would be valid for movement to the surface + exceeds this value the simulation will be interrupted + +* **toleranceMoveSurfaceReactionToCore** + + A surface reaction will be pulled into the bulk core if its dynamics number + ever exceeds this value. Note this is done on the fly during simulation. + +* **toleranceMoveSurfaceSpeciesToCore** + + A surface species will be pulled into the bulk core if it's rate ratio + ever exceeds this value. Note this is done on the fly during simulation. diff --git a/documentation/source/theory/rmg/index.rst b/documentation/source/theory/rmg/index.rst index ba074c472b..92d1d66ec0 100644 --- a/documentation/source/theory/rmg/index.rst +++ b/documentation/source/theory/rmg/index.rst @@ -10,6 +10,7 @@ RMG Theory Guide enlarger prune + dynamics * :ref:`genindex` * :ref:`modindex` diff --git a/documentation/source/theory/rmg/surfacediagram.png b/documentation/source/theory/rmg/surfacediagram.png new file mode 100644 index 0000000000..e71ea5f6c7 Binary files /dev/null and b/documentation/source/theory/rmg/surfacediagram.png differ diff --git a/documentation/source/users/cantherm/input.rst b/documentation/source/users/cantherm/input.rst index a7c6976e70..563ca17f4e 100644 --- a/documentation/source/users/cantherm/input.rst +++ b/documentation/source/users/cantherm/input.rst @@ -76,6 +76,7 @@ Model Chemistry AEC BC SOC ``'BMK/cbsb7'`` v v v ``'BMK/6-311G(2d,d,p)'`` v v v ``'B3LYP/6-311+G(3df,2p)'`` v +``'B3LYP/6-31G**'`` v v ================================================ ===== ==== ==== Notes: @@ -122,8 +123,8 @@ Parameter Required? Description ``externalSymmetry`` yes The external symmetry number for rotation ``spinMultiplicity`` yes The ground-state spin multiplicity (degeneracy) ``opticalIsomers`` yes The number of optical isomers of the species -``energy`` yes The ground-state 0 K atomization energy in Hartree (without zero-point energy) - **or** +``energy`` yes The ground-state 0 K atomization energy in Hartree + (without zero-point energy) **or** The path to the quantum chemistry output file containing the energy ``geometry`` yes The path to the quantum chemistry output file containing the optimized geometry ``frequencies`` yes The path to the quantum chemistry output file containing the computed frequencies @@ -170,7 +171,7 @@ they can specify the path to a quantum chemistry calculation output file that co In this example, the ``CBS-QB3`` energy is obtained from a Gaussian log file, while the ``Klip_2`` energy is specified directly. The energy used will depend on what ``modelChemistry()`` was specified in the input file. CanTherm can parse the energy from -a ``GaussianLog``, ``MoleProLog`` or ``QchemLog``. +a ``GaussianLog``, ``MolproLog`` or ``QchemLog``. The input to the remaining parameters, ``geometry``, ``frequencies`` and ``rotors``, will depend on if hindered/free rotors are included. Both cases are described below. @@ -266,6 +267,24 @@ As noted above, ``scanLog`` can either point to a ``GaussianLog``, ``QchemLog`` The ``Energy`` can be in units of ``kJ/mol``, ``J/mol``, ``cal/mol``, ``kcal/mol``, ``cm^-1`` or ``hartree``. +The ``symmetry`` parameter will usually equal either 1, 2 or 3. Below are examples of internal rotor scans with these commonly encountered symmetry numbers. First, ``symmetry = 3``: + +.. image:: symmetry_3_example.png + +Internal rotation of a methyl group is a common example of a hindered rotor with ``symmetry = 3``, such as the one above. As shown, all three minima (and maxima) have identical energies, hence ``symmetry = 3``. + +Similarly, if there are only two minima along the internal rotor scan, and both have identical energy, then ``symmetry = 2``, as in the example below: + +.. image:: symmetry_2_example.png + +If any of the energy minima in an internal rotor scan are not identical, then the rotor has no symmetry (``symmetry = 1``), as in the example below: + +.. image:: symmetry_1_example.png + +For the example above there are 3 local energy minima, 2 of which are identical to each other. However, the 3rd minima is different from the other 2, therefore this internal rotor has no symmetry. + +For practical purposes, when determining the symmetry number for a given hindered rotor simply check if the internal rotor scan looks like the ``symmetry = 2`` or ``3`` examples above. If it doesn’t, then most likely ``symmetry = 1``. + Each :class:`FreeRotor()` object requires the following parameters: ====================== ========================================================= diff --git a/documentation/source/users/cantherm/input_pdep.rst b/documentation/source/users/cantherm/input_pdep.rst index 87ba226703..30069164db 100644 --- a/documentation/source/users/cantherm/input_pdep.rst +++ b/documentation/source/users/cantherm/input_pdep.rst @@ -170,8 +170,8 @@ Parameter Required? Description ``externalSymmetry`` yes The external symmetry number for rotation ``spinMultiplicity`` yes The ground-state spin multiplicity (degeneracy) ``opticalIsomers`` yes The number of optical isomers of the species -``energy`` yes The ground-state 0 K atomization energy in Hartree (without zero-point energy) - **or** +``energy`` yes The ground-state 0 K atomization energy in Hartree + (without zero-point energy) **or** The path to the quantum chemistry output file containing the energy ``geometry`` yes The path to the quantum chemistry output file containing the optimized geometry ``frequencies`` yes The path to the quantum chemistry output file containing the computed frequencies @@ -218,7 +218,7 @@ they can specify the path to a quantum chemistry calculation output file that co In this example, the ``CBS-QB3`` energy is obtained from a Gaussian log file, while the ``Klip_2`` energy is specified directly. The energy used will depend on what ``modelChemistry()`` was specified in the input file. CanTherm can parse the energy from -a ``GaussianLog``, ``MoleProLog`` or ``QchemLog``. +a ``GaussianLog``, ``MolproLog`` or ``QchemLog``. The input to the remaining parameters, ``geometry``, ``frequencies`` and ``rotors``, will depend on if hindered/free rotors are included. Both cases are described below. @@ -314,6 +314,24 @@ As noted above, ``scanLog`` can either point to a ``GaussianLog``, ``QchemLog`` The ``Energy`` can be in units of ``kJ/mol``, ``J/mol``, ``cal/mol``, ``kcal/mol``, ``cm^-1`` or ``hartree``. +The ``symmetry`` parameter will usually equal either 1, 2 or 3. Below are examples of internal rotor scans with these commonly encountered symmetry numbers. First, ``symmetry = 3``: + +.. image:: symmetry_3_example.png + +Internal rotation of a methyl group is a common example of a hindered rotor with ``symmetry = 3``, such as the one above. As shown, all three minima (and maxima) have identical energies, hence ``symmetry = 3``. + +Similarly, if there are only two minima along the internal rotor scan, and both have identical energy, then ``symmetry = 2``, as in the example below: + +.. image:: symmetry_2_example.png + +If any of the energy minima in an internal rotor scan are not identical, then the rotor has no symmetry (``symmetry = 1``), as in the example below: + +.. image:: symmetry_1_example.png + +For the example above there are 3 local energy minima, 2 of which are identical to each other. However, the 3rd minima is different from the other 2, therefore this internal rotor has no symmetry. + +For practical purposes, when determining the symmetry number for a given hindered rotor simply check if the internal rotor scan looks like the ``symmetry = 2`` or ``3`` examples above. If it doesn’t, then most likely ``symmetry = 1``. + Each :class:`FreeRotor()` object requires the following parameters: ====================== ========================================================= diff --git a/documentation/source/users/cantherm/output.rst b/documentation/source/users/cantherm/output.rst index a3f585c9b1..7876cf273e 100644 --- a/documentation/source/users/cantherm/output.rst +++ b/documentation/source/users/cantherm/output.rst @@ -57,3 +57,11 @@ something unexpected has occurred. The ``examples/cantherm`` directory contains both CanTherm input files and the resulting output files. + +Species Dictionary +================== + +Any species that had the ``thermo()`` method called and had the structure defined in the cantherm +input file will also have an RMG style adjacency list representation in ``species_dictionary.txt``. +This allows the user to input the corresponding thermo and kinetics into RMG in various ways +described in the RMG user guide. \ No newline at end of file diff --git a/documentation/source/users/cantherm/symmetry_1_example.png b/documentation/source/users/cantherm/symmetry_1_example.png new file mode 100644 index 0000000000..0d60afba78 Binary files /dev/null and b/documentation/source/users/cantherm/symmetry_1_example.png differ diff --git a/documentation/source/users/cantherm/symmetry_2_example.png b/documentation/source/users/cantherm/symmetry_2_example.png new file mode 100644 index 0000000000..61521c981e Binary files /dev/null and b/documentation/source/users/cantherm/symmetry_2_example.png differ diff --git a/documentation/source/users/cantherm/symmetry_3_example.png b/documentation/source/users/cantherm/symmetry_3_example.png new file mode 100644 index 0000000000..9cf5a66eba Binary files /dev/null and b/documentation/source/users/cantherm/symmetry_3_example.png differ diff --git a/documentation/source/users/rmg/credits.rst b/documentation/source/users/rmg/credits.rst index 8d75a1f99b..ba9d0d4563 100755 --- a/documentation/source/users/rmg/credits.rst +++ b/documentation/source/users/rmg/credits.rst @@ -16,9 +16,9 @@ Current Developers: (rmg_dev@mit.edu) - Dr. Alon G. Dana - Mark Goldman - Kehang Han -- Max Liu -- Belinda Slakman -- Nathan Yee +- Matt Johnson +- Mengjie Liu +- Mark Payne Previous Developers: @@ -31,9 +31,12 @@ Previous Developers: - Dr. Fariba S. Khanshan - Victor Lambert - Dr. Shamel S. Merchant +- Dr. Belinda Slakman - Sean Troiano - Dr. Aaron Vandeputte - Dr. Nick M. Vandewiele +- Dr. Nathan Yee +- Peng Zhang *********** diff --git a/documentation/source/users/rmg/database/images/kinetics_families/1+2_Cycloaddition.png b/documentation/source/users/rmg/database/images/kinetics_families/1+2_Cycloaddition.png index c57976e89c..cc6bbe24d9 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/1+2_Cycloaddition.png and b/documentation/source/users/rmg/database/images/kinetics_families/1+2_Cycloaddition.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/1,2-Birad_to_alkene.png b/documentation/source/users/rmg/database/images/kinetics_families/1,2-Birad_to_alkene.png index 9b88365b5a..ea17be3e4a 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/1,2-Birad_to_alkene.png and b/documentation/source/users/rmg/database/images/kinetics_families/1,2-Birad_to_alkene.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/1,2_Insertion_CO.png b/documentation/source/users/rmg/database/images/kinetics_families/1,2_Insertion_CO.png index 3fa7d99ab6..97e9662d6b 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/1,2_Insertion_CO.png and b/documentation/source/users/rmg/database/images/kinetics_families/1,2_Insertion_CO.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/1,2_Insertion_carbene.png b/documentation/source/users/rmg/database/images/kinetics_families/1,2_Insertion_carbene.png index 2a8eb03761..82b050390c 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/1,2_Insertion_carbene.png and b/documentation/source/users/rmg/database/images/kinetics_families/1,2_Insertion_carbene.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/1,2_shiftS.png b/documentation/source/users/rmg/database/images/kinetics_families/1,2_shiftS.png index 289d224fff..e6d5817da7 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/1,2_shiftS.png and b/documentation/source/users/rmg/database/images/kinetics_families/1,2_shiftS.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/1,3_Insertion_CO2.png b/documentation/source/users/rmg/database/images/kinetics_families/1,3_Insertion_CO2.png index b8ca2d4c6e..90af820ab3 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/1,3_Insertion_CO2.png and b/documentation/source/users/rmg/database/images/kinetics_families/1,3_Insertion_CO2.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/1,3_Insertion_ROR.png b/documentation/source/users/rmg/database/images/kinetics_families/1,3_Insertion_ROR.png index 2c96479fe0..2b4cf4a2cb 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/1,3_Insertion_ROR.png and b/documentation/source/users/rmg/database/images/kinetics_families/1,3_Insertion_ROR.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/1,3_Insertion_RSR.png b/documentation/source/users/rmg/database/images/kinetics_families/1,3_Insertion_RSR.png index 652373bab3..e30f9af13c 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/1,3_Insertion_RSR.png and b/documentation/source/users/rmg/database/images/kinetics_families/1,3_Insertion_RSR.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/1,4_Cyclic_birad_scission.png b/documentation/source/users/rmg/database/images/kinetics_families/1,4_Cyclic_birad_scission.png index cd8528e075..5d6d107c0d 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/1,4_Cyclic_birad_scission.png and b/documentation/source/users/rmg/database/images/kinetics_families/1,4_Cyclic_birad_scission.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/1,4_Linear_birad_scission.png b/documentation/source/users/rmg/database/images/kinetics_families/1,4_Linear_birad_scission.png index f518691f83..0053543011 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/1,4_Linear_birad_scission.png and b/documentation/source/users/rmg/database/images/kinetics_families/1,4_Linear_birad_scission.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/2+2_cycloaddition_CCO.png b/documentation/source/users/rmg/database/images/kinetics_families/2+2_cycloaddition_CCO.png index b7262f078d..7db1159ca6 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/2+2_cycloaddition_CCO.png and b/documentation/source/users/rmg/database/images/kinetics_families/2+2_cycloaddition_CCO.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/2+2_cycloaddition_CO.png b/documentation/source/users/rmg/database/images/kinetics_families/2+2_cycloaddition_CO.png index aac04d0a2c..5bbe3bd0c7 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/2+2_cycloaddition_CO.png and b/documentation/source/users/rmg/database/images/kinetics_families/2+2_cycloaddition_CO.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/2+2_cycloaddition_CS.png b/documentation/source/users/rmg/database/images/kinetics_families/2+2_cycloaddition_CS.png new file mode 100644 index 0000000000..d83a675662 Binary files /dev/null and b/documentation/source/users/rmg/database/images/kinetics_families/2+2_cycloaddition_CS.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/2+2_cycloaddition_Cd.png b/documentation/source/users/rmg/database/images/kinetics_families/2+2_cycloaddition_Cd.png index d19227beb4..52f4b62258 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/2+2_cycloaddition_Cd.png and b/documentation/source/users/rmg/database/images/kinetics_families/2+2_cycloaddition_Cd.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/6_membered_central_C-C_shift.png b/documentation/source/users/rmg/database/images/kinetics_families/6_membered_central_C-C_shift.png new file mode 100644 index 0000000000..947f02b362 Binary files /dev/null and b/documentation/source/users/rmg/database/images/kinetics_families/6_membered_central_C-C_shift.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/Birad_R_Recombination.png b/documentation/source/users/rmg/database/images/kinetics_families/Birad_R_Recombination.png new file mode 100644 index 0000000000..0adc4069f3 Binary files /dev/null and b/documentation/source/users/rmg/database/images/kinetics_families/Birad_R_Recombination.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/Birad_recombination.png b/documentation/source/users/rmg/database/images/kinetics_families/Birad_recombination.png index 84f5e15907..124dcfee25 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/Birad_recombination.png and b/documentation/source/users/rmg/database/images/kinetics_families/Birad_recombination.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/CO_Disproportionation.png b/documentation/source/users/rmg/database/images/kinetics_families/CO_Disproportionation.png new file mode 100644 index 0000000000..1a6d984405 Binary files /dev/null and b/documentation/source/users/rmg/database/images/kinetics_families/CO_Disproportionation.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/Concerted_Intra_Diels_alder_monocyclic_1,2_shiftH.png b/documentation/source/users/rmg/database/images/kinetics_families/Concerted_Intra_Diels_alder_monocyclic_1,2_shiftH.png new file mode 100644 index 0000000000..a8732830bb Binary files /dev/null and b/documentation/source/users/rmg/database/images/kinetics_families/Concerted_Intra_Diels_alder_monocyclic_1,2_shiftH.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/Cyclic_Ether_Formation.png b/documentation/source/users/rmg/database/images/kinetics_families/Cyclic_Ether_Formation.png index 430cac700a..e7189421f1 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/Cyclic_Ether_Formation.png and b/documentation/source/users/rmg/database/images/kinetics_families/Cyclic_Ether_Formation.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/Cyclic_Thioether_Formation.png b/documentation/source/users/rmg/database/images/kinetics_families/Cyclic_Thioether_Formation.png new file mode 100644 index 0000000000..9738b9bb85 Binary files /dev/null and b/documentation/source/users/rmg/database/images/kinetics_families/Cyclic_Thioether_Formation.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/Cyclopentadiene_scission.png b/documentation/source/users/rmg/database/images/kinetics_families/Cyclopentadiene_scission.png new file mode 100644 index 0000000000..aa4c65c53a Binary files /dev/null and b/documentation/source/users/rmg/database/images/kinetics_families/Cyclopentadiene_scission.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/Diels_alder_addition.png b/documentation/source/users/rmg/database/images/kinetics_families/Diels_alder_addition.png index 1f03a7d857..64ace091c3 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/Diels_alder_addition.png and b/documentation/source/users/rmg/database/images/kinetics_families/Diels_alder_addition.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/Disproportionation.png b/documentation/source/users/rmg/database/images/kinetics_families/Disproportionation.png index a98efeea92..69e8b26842 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/Disproportionation.png and b/documentation/source/users/rmg/database/images/kinetics_families/Disproportionation.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/HO2_Elimination_from_PeroxyRadical.png b/documentation/source/users/rmg/database/images/kinetics_families/HO2_Elimination_from_PeroxyRadical.png index 7d0c38e8d2..2b20bdcbcb 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/HO2_Elimination_from_PeroxyRadical.png and b/documentation/source/users/rmg/database/images/kinetics_families/HO2_Elimination_from_PeroxyRadical.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/H_Abstraction.png b/documentation/source/users/rmg/database/images/kinetics_families/H_Abstraction.png index 06580069a5..f09ad770f2 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/H_Abstraction.png and b/documentation/source/users/rmg/database/images/kinetics_families/H_Abstraction.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/H_shift_cyclopentadiene.png b/documentation/source/users/rmg/database/images/kinetics_families/H_shift_cyclopentadiene.png deleted file mode 100644 index 80aa0f7b34..0000000000 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/H_shift_cyclopentadiene.png and /dev/null differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/Intra_2+2_cycloaddition_Cd.png b/documentation/source/users/rmg/database/images/kinetics_families/Intra_2+2_cycloaddition_Cd.png new file mode 100644 index 0000000000..b2b3289aca Binary files /dev/null and b/documentation/source/users/rmg/database/images/kinetics_families/Intra_2+2_cycloaddition_Cd.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/Intra_5_membered_conjugated_C=C_C=C_addition.png b/documentation/source/users/rmg/database/images/kinetics_families/Intra_5_membered_conjugated_C=C_C=C_addition.png new file mode 100644 index 0000000000..54f4894318 Binary files /dev/null and b/documentation/source/users/rmg/database/images/kinetics_families/Intra_5_membered_conjugated_C=C_C=C_addition.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/Intra_Diels_alder.png b/documentation/source/users/rmg/database/images/kinetics_families/Intra_Diels_alder.png deleted file mode 100644 index 9f5d45a74f..0000000000 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/Intra_Diels_alder.png and /dev/null differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/Intra_Diels_alder_monocyclic.png b/documentation/source/users/rmg/database/images/kinetics_families/Intra_Diels_alder_monocyclic.png new file mode 100644 index 0000000000..e469b31a77 Binary files /dev/null and b/documentation/source/users/rmg/database/images/kinetics_families/Intra_Diels_alder_monocyclic.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/Intra_Disproportionation.png b/documentation/source/users/rmg/database/images/kinetics_families/Intra_Disproportionation.png index d6d50012ea..6def6795ca 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/Intra_Disproportionation.png and b/documentation/source/users/rmg/database/images/kinetics_families/Intra_Disproportionation.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/Intra_RH_Add_Endocyclic.png b/documentation/source/users/rmg/database/images/kinetics_families/Intra_RH_Add_Endocyclic.png index 18956f3f5f..d292851e99 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/Intra_RH_Add_Endocyclic.png and b/documentation/source/users/rmg/database/images/kinetics_families/Intra_RH_Add_Endocyclic.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/Intra_RH_Add_Exocyclic.png b/documentation/source/users/rmg/database/images/kinetics_families/Intra_RH_Add_Exocyclic.png index 25bb064c7a..b26964707d 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/Intra_RH_Add_Exocyclic.png and b/documentation/source/users/rmg/database/images/kinetics_families/Intra_RH_Add_Exocyclic.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/Intra_R_Add_Endocyclic.png b/documentation/source/users/rmg/database/images/kinetics_families/Intra_R_Add_Endocyclic.png index 7fad539dc4..183898f825 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/Intra_R_Add_Endocyclic.png and b/documentation/source/users/rmg/database/images/kinetics_families/Intra_R_Add_Endocyclic.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/Intra_R_Add_ExoTetCyclic.png b/documentation/source/users/rmg/database/images/kinetics_families/Intra_R_Add_ExoTetCyclic.png index 2bd4094d0c..e037b9f149 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/Intra_R_Add_ExoTetCyclic.png and b/documentation/source/users/rmg/database/images/kinetics_families/Intra_R_Add_ExoTetCyclic.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/Intra_R_Add_Exocyclic.png b/documentation/source/users/rmg/database/images/kinetics_families/Intra_R_Add_Exocyclic.png index d07c395035..582867e1af 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/Intra_R_Add_Exocyclic.png and b/documentation/source/users/rmg/database/images/kinetics_families/Intra_R_Add_Exocyclic.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/Intra_Retro_Diels_alder_bicyclic.png b/documentation/source/users/rmg/database/images/kinetics_families/Intra_Retro_Diels_alder_bicyclic.png new file mode 100644 index 0000000000..8e0c246ea3 Binary files /dev/null and b/documentation/source/users/rmg/database/images/kinetics_families/Intra_Retro_Diels_alder_bicyclic.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/Intra_ene_reaction.png b/documentation/source/users/rmg/database/images/kinetics_families/Intra_ene_reaction.png new file mode 100644 index 0000000000..2ddf95fbf8 Binary files /dev/null and b/documentation/source/users/rmg/database/images/kinetics_families/Intra_ene_reaction.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/Korcek_step1.png b/documentation/source/users/rmg/database/images/kinetics_families/Korcek_step1.png index 0b3d0edc13..a9a476ec1c 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/Korcek_step1.png and b/documentation/source/users/rmg/database/images/kinetics_families/Korcek_step1.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/Korcek_step2.png b/documentation/source/users/rmg/database/images/kinetics_families/Korcek_step2.png index a6acd44272..b56e9d2166 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/Korcek_step2.png and b/documentation/source/users/rmg/database/images/kinetics_families/Korcek_step2.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/Oa_R_Recombination.png b/documentation/source/users/rmg/database/images/kinetics_families/Oa_R_Recombination.png deleted file mode 100644 index db1aae1252..0000000000 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/Oa_R_Recombination.png and /dev/null differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/R_Addition_COm.png b/documentation/source/users/rmg/database/images/kinetics_families/R_Addition_COm.png index 7e11303340..7b5427475d 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/R_Addition_COm.png and b/documentation/source/users/rmg/database/images/kinetics_families/R_Addition_COm.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/R_Addition_CSm.png b/documentation/source/users/rmg/database/images/kinetics_families/R_Addition_CSm.png index b20c71bc02..9522d72eac 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/R_Addition_CSm.png and b/documentation/source/users/rmg/database/images/kinetics_families/R_Addition_CSm.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/R_Addition_MultipleBond.png b/documentation/source/users/rmg/database/images/kinetics_families/R_Addition_MultipleBond.png index b8d1f1e873..adc9416930 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/R_Addition_MultipleBond.png and b/documentation/source/users/rmg/database/images/kinetics_families/R_Addition_MultipleBond.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/R_Recombination.png b/documentation/source/users/rmg/database/images/kinetics_families/R_Recombination.png index ae232ec737..6db7b5397f 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/R_Recombination.png and b/documentation/source/users/rmg/database/images/kinetics_families/R_Recombination.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/Singlet_Carbene_Intra_Disproportionation.png b/documentation/source/users/rmg/database/images/kinetics_families/Singlet_Carbene_Intra_Disproportionation.png new file mode 100644 index 0000000000..4bed86959b Binary files /dev/null and b/documentation/source/users/rmg/database/images/kinetics_families/Singlet_Carbene_Intra_Disproportionation.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/Singlet_Val6_to_triplet.png b/documentation/source/users/rmg/database/images/kinetics_families/Singlet_Val6_to_triplet.png new file mode 100644 index 0000000000..add6bc8788 Binary files /dev/null and b/documentation/source/users/rmg/database/images/kinetics_families/Singlet_Val6_to_triplet.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/SubstitutionS.png b/documentation/source/users/rmg/database/images/kinetics_families/SubstitutionS.png index b94449df1e..578108478c 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/SubstitutionS.png and b/documentation/source/users/rmg/database/images/kinetics_families/SubstitutionS.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/Substitution_O.png b/documentation/source/users/rmg/database/images/kinetics_families/Substitution_O.png index 4298165b85..152a811e57 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/Substitution_O.png and b/documentation/source/users/rmg/database/images/kinetics_families/Substitution_O.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/intra_H_migration.png b/documentation/source/users/rmg/database/images/kinetics_families/intra_H_migration.png index 71634f9c10..5dd5cbb42a 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/intra_H_migration.png and b/documentation/source/users/rmg/database/images/kinetics_families/intra_H_migration.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/intra_NO2_ONO_conversion.png b/documentation/source/users/rmg/database/images/kinetics_families/intra_NO2_ONO_conversion.png index 8023fd16e9..a0ac930c4d 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/intra_NO2_ONO_conversion.png and b/documentation/source/users/rmg/database/images/kinetics_families/intra_NO2_ONO_conversion.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/intra_OH_migration.png b/documentation/source/users/rmg/database/images/kinetics_families/intra_OH_migration.png index 6d43a5b1eb..7e373bfa44 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/intra_OH_migration.png and b/documentation/source/users/rmg/database/images/kinetics_families/intra_OH_migration.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/intra_substitutionCS_cyclization.png b/documentation/source/users/rmg/database/images/kinetics_families/intra_substitutionCS_cyclization.png index 4258929237..d10cd056e6 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/intra_substitutionCS_cyclization.png and b/documentation/source/users/rmg/database/images/kinetics_families/intra_substitutionCS_cyclization.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/intra_substitutionCS_isomerization.png b/documentation/source/users/rmg/database/images/kinetics_families/intra_substitutionCS_isomerization.png index d098cd0de2..a4644cc7eb 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/intra_substitutionCS_isomerization.png and b/documentation/source/users/rmg/database/images/kinetics_families/intra_substitutionCS_isomerization.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/intra_substitutionS_cyclization.png b/documentation/source/users/rmg/database/images/kinetics_families/intra_substitutionS_cyclization.png index 2bb887d16a..9eb5e3b6e8 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/intra_substitutionS_cyclization.png and b/documentation/source/users/rmg/database/images/kinetics_families/intra_substitutionS_cyclization.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/intra_substitutionS_isomerization.png b/documentation/source/users/rmg/database/images/kinetics_families/intra_substitutionS_isomerization.png index 8e7db6ebe4..297b5380ad 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/intra_substitutionS_isomerization.png and b/documentation/source/users/rmg/database/images/kinetics_families/intra_substitutionS_isomerization.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/ketoenol.png b/documentation/source/users/rmg/database/images/kinetics_families/ketoenol.png index 2c8147b06a..f0034fcba3 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/ketoenol.png and b/documentation/source/users/rmg/database/images/kinetics_families/ketoenol.png differ diff --git a/documentation/source/users/rmg/database/images/kinetics_families/lone_electron_pair_bond.png b/documentation/source/users/rmg/database/images/kinetics_families/lone_electron_pair_bond.png index 3c8c4e914f..d2a2d06408 100644 Binary files a/documentation/source/users/rmg/database/images/kinetics_families/lone_electron_pair_bond.png and b/documentation/source/users/rmg/database/images/kinetics_families/lone_electron_pair_bond.png differ diff --git a/documentation/source/users/rmg/database/introduction.rst b/documentation/source/users/rmg/database/introduction.rst index 1e5ab197dd..43b1add114 100644 --- a/documentation/source/users/rmg/database/introduction.rst +++ b/documentation/source/users/rmg/database/introduction.rst @@ -155,7 +155,7 @@ table below shows all atoms types in RMG. +----------+-------------------+------------------------------------------------------------------------------------------------------------------+ |S4s |Sulfur |One lone pair, up to four single bonds | +----------+-------------------+------------------------------------------------------------------------------------------------------------------+ -|S4sc |Sulfur |One lone pair, up to five single bonds, charged -1/+1 | +|S4sc |Sulfur |One lone pair, up to five single bonds, charged -1/+1 | +----------+-------------------+------------------------------------------------------------------------------------------------------------------+ |S4d |Sulfur |One lone pair, one double bond, up to two single bonds | +----------+-------------------+------------------------------------------------------------------------------------------------------------------+ @@ -189,6 +189,14 @@ table below shows all atoms types in RMG. +----------+-------------------+------------------------------------------------------------------------------------------------------------------+ |S6tdc |Sulfur |No lone pairs, one to two triple bonds, up to two double bonds, up to four single bonds, charged -1/-1 | +----------+-------------------+------------------------------------------------------------------------------------------------------------------+ +|Cl |Chlorine |No requirements | ++----------+-------------------+------------------------------------------------------------------------------------------------------------------+ +|Cl1s |Chlorine |Three lone pairs, zero to one single bonds | ++----------+-------------------+------------------------------------------------------------------------------------------------------------------+ +|I |Iodine |No requirements | ++----------+-------------------+------------------------------------------------------------------------------------------------------------------+ +|I1s |Iodine |Three lone pairs, zero to one single bonds | ++----------+-------------------+------------------------------------------------------------------------------------------------------------------+ |He |Helium |No requirements, nonreactive | +----------+-------------------+------------------------------------------------------------------------------------------------------------------+ |Ne |Neon |No requirements, nonreactive | @@ -214,12 +222,12 @@ or products that are forbidden. These groups are stored in in the file located To ban certain specific pathways in the kinetics families, a `forbidden` group must be created, like the following group -in the ``intra_H_migration`` family :: +in the ``intra_H_migration`` family: :: forbidden( label = "bridged56_1254", - group = - """"""" + group = + """ 1 *1 C 1 {2,S} {6,S} 2 *4 C 0 {1,S} {3,S} {7,S} 3 C 0 {2,S} {4,S} @@ -230,9 +238,9 @@ in the ``intra_H_migration`` family :: 8 *3 H 0 {4,S} """, shortDesc = u"""""", - longDesc = + longDesc = u""" - + """, ) @@ -241,6 +249,42 @@ specific kinetics family's folder ``RMG-database/input/kinetics/family_name/`` alongside normal group entries. The starred atoms in the forbidden group ban the specified reaction recipe from occurring in matched products and reactants. +In addition for forbidding groups, there is the option of forbidding specific +molecules or species. Forbidding a molecule will prevent that exact structure +from being generated, while forbidding a species will prevent any of its resonance +structures from being generated. To specify a forbidden molecule or species, simply +replace the ``group`` keyword with ``molecule`` or ``species``: :: + + # This forbids a molecule + forbidden( + label = "C_quintet", + molecule = + """ + multiplicity 5 + 1 C u4 p0 + """, + shortDesc = u"""""", + longDesc = + u""" + + """, + ) + + # This forbids a species + forbidden( + label = "C_quintet", + species = + """ + multiplicity 5 + 1 C u4 p0 + """, + shortDesc = u"""""", + longDesc = + u""" + + """, + ) + Hierarchical Trees ------------------ Groups are ordered into the nodes of a hierarchical trees which is written diff --git a/documentation/source/users/rmg/database/kinetics.rst b/documentation/source/users/rmg/database/kinetics.rst index 41c37b977d..f9cf0a6d9b 100644 --- a/documentation/source/users/rmg/database/kinetics.rst +++ b/documentation/source/users/rmg/database/kinetics.rst @@ -211,234 +211,182 @@ Each reaction family contains the files: * training.py containing a training set for the family * rules.py containing kinetic parameters for rules -There are currently 45 reaction families in RMG: - -**1+2_Cycloaddition** - -.. image:: images/kinetics_families/1+2_Cycloaddition.png - :scale: 40% - -**1,2-Birad_to_alkene** - -.. image:: images/kinetics_families/1,2-Birad_to_alkene.png - :scale: 40% - -**1,2_Insertion_carbene** - -.. image:: images/kinetics_families/1,2_Insertion_carbene.png - :scale: 40% - -**1,2_Insertion_CO** - -.. image:: images/kinetics_families/1,2_Insertion_CO.png - :scale: 40% - -**1,2_shiftS** - -.. image:: images/kinetics_families/1,2_shiftS.png - :scale: 40% - -**1,3_Insertion_CO2** - -.. image:: images/kinetics_families/1,3_Insertion_CO2.png - :scale: 40% - -**1,3_Insertion_ROR** - -.. image:: images/kinetics_families/1,3_Insertion_ROR.png - :scale: 40% - -**1,3_Insertion_RSR** - -.. image:: images/kinetics_families/1,3_Insertion_RSR.png - :scale: 40% - -**1,4_Cyclic_birad_scission** - -.. image:: images/kinetics_families/1,4_Cyclic_birad_scission.png - :scale: 40% - -**1,4_Linear_birad_scission** - -.. image:: images/kinetics_families/1,4_Linear_birad_scission.png - :scale: 40% - -**2+2_cycloaddition_CCO** - -.. image:: images/kinetics_families/2+2_cycloaddition_CCO.png - :scale: 40% - -**2+2_cycloaddition_Cd** - -.. image:: images/kinetics_families/2+2_cycloaddition_Cd.png - :scale: 40% - -**2+2_cycloaddition_CO** - -.. image:: images/kinetics_families/2+2_cycloaddition_CO.png - :scale: 40% - -**Birad_recombination** - -.. image:: images/kinetics_families/Birad_recombination.png - :scale: 40% - -**Cyclic_Ether_Formation** - -.. image:: images/kinetics_families/Cyclic_Ether_Formation.png - :scale: 40% - -**Diels_alder_addition** - -.. image:: images/kinetics_families/Diels_alder_addition.png - :scale: 40% - -**Disproportionation** - -.. image:: images/kinetics_families/Disproportionation.png - :scale: 40% - -**H_Abstraction** - -.. image:: images/kinetics_families/H_Abstraction.png - :scale: 40% - -**H_shift_cyclopentadiene** - -.. image:: images/kinetics_families/H_shift_cyclopentadiene.png - :scale: 40% - -**HO2_Elimination_from_PeroxyRadical** - -.. image:: images/kinetics_families/HO2_Elimination_from_PeroxyRadical.png - :scale: 40% - -**Intra_Diels_alder** - -.. image:: images/kinetics_families/Intra_Diels_alder.png - :scale: 40% - -**Intra_Disproportionation** - -.. image:: images/kinetics_families/Intra_Disproportionation.png - :scale: 40% - -**intra_H_migration** - -.. image:: images/kinetics_families/intra_H_migration.png - :scale: 40% - -**intra_NO2_ONO_conversion** - -.. image:: images/kinetics_families/intra_NO2_ONO_conversion.png - :scale: 40% - -**intra_OH_migration** - -.. image:: images/kinetics_families/intra_OH_migration.png - :scale: 40% - -**Intra_R_Add_Endocyclic** - -.. image:: images/kinetics_families/Intra_R_Add_Endocyclic.png - :scale: 40% - -**Intra_R_Add_Exocyclic** - -.. image:: images/kinetics_families/Intra_R_Add_Exocyclic.png - :scale: 40% - -**Intra_R_Add_ExoTetCyclic** - -.. image:: images/kinetics_families/Intra_R_Add_ExoTetCyclic.png - :scale: 40% - -**Intra_RH_Add_Endocyclic** - -.. image:: images/kinetics_families/Intra_RH_Add_Endocyclic.png - :scale: 40% - -**Intra_RH_Add_Exocyclic** - -.. image:: images/kinetics_families/Intra_RH_Add_Exocyclic.png - :scale: 40% - -**intra_substitutionCS_cyclization** - -.. image:: images/kinetics_families/intra_substitutionCS_cyclization.png - :scale: 40% - -**intra_substitutionCS_isomerization** - -.. image:: images/kinetics_families/intra_substitutionCS_isomerization.png - :scale: 40% - -**intra_substitutionS_cyclization** - -.. image:: images/kinetics_families/intra_substitutionS_cyclization.png - :scale: 40% - -**intra_substitutionS_isomerization** - -.. image:: images/kinetics_families/intra_substitutionS_isomerization.png - :scale: 40% - -**ketoenol** - -.. image:: images/kinetics_families/ketoenol.png - :scale: 40% - -**Korcek_step1** - -.. image:: images/kinetics_families/Korcek_step1.png - :scale: 40% - -**Korcek_step2** - -.. image:: images/kinetics_families/Korcek_step2.png - :scale: 40% - -**lone_electron_pair_bond** - -.. image:: images/kinetics_families/lone_electron_pair_bond.png - :scale: 40% - -**Oa_R_Recombination** - -.. image:: images/kinetics_families/Oa_R_Recombination.png - :scale: 40% - -**R_Addition_COm** - -.. image:: images/kinetics_families/R_Addition_COm.png - :scale: 40% - -**R_Addition_CSm** - -.. image:: images/kinetics_families/R_Addition_CSm.png - :scale: 40% - -**R_Addition_MultipleBond** - -.. image:: images/kinetics_families/R_Addition_MultipleBond.png - :scale: 40% - -**R_Recombination** - -.. image:: images/kinetics_families/R_Recombination.png - :scale: 40% - -**Substitution_O** - -.. image:: images/kinetics_families/Substitution_O.png - :scale: 40% - -**SubstitutionS** - -.. image:: images/kinetics_families/SubstitutionS.png - :scale: 40% - - +There are currently 58 reaction families in RMG: + +.. |f01| image:: images/kinetics_families/1,2-Birad_to_alkene.png + :scale: 25% +.. |f02| image:: images/kinetics_families/1,2_Insertion_carbene.png + :scale: 25% +.. |f03| image:: images/kinetics_families/1,2_Insertion_CO.png + :scale: 25% +.. |f05| image:: images/kinetics_families/1,2_shiftS.png + :scale: 25% +.. |f06| image:: images/kinetics_families/1,3_Insertion_CO2.png + :scale: 25% +.. |f07| image:: images/kinetics_families/1,3_Insertion_ROR.png + :scale: 25% +.. |f08| image:: images/kinetics_families/1,3_Insertion_RSR.png + :scale: 25% +.. |f09| image:: images/kinetics_families/1,4_Cyclic_birad_scission.png + :scale: 25% +.. |f10| image:: images/kinetics_families/1,4_Linear_birad_scission.png + :scale: 25% +.. |f11| image:: images/kinetics_families/1+2_Cycloaddition.png + :scale: 25% +.. |f12| image:: images/kinetics_families/2+2_cycloaddition_CCO.png + :scale: 25% +.. |f13| image:: images/kinetics_families/2+2_cycloaddition_Cd.png + :scale: 25% +.. |f14| image:: images/kinetics_families/2+2_cycloaddition_CO.png + :scale: 25% +.. |f15| image:: images/kinetics_families/2+2_cycloaddition_CS.png + :scale: 25% +.. |f16| image:: images/kinetics_families/6_membered_central_C-C_shift.png + :scale: 25% +.. |f17| image:: images/kinetics_families/Birad_recombination.png + :scale: 25% +.. |f18| image:: images/kinetics_families/Birad_R_Recombination.png + :scale: 25% +.. |f19| image:: images/kinetics_families/CO_Disproportionation.png + :scale: 25% +.. |f20| image:: images/kinetics_families/Concerted_Intra_Diels_alder_monocyclic_1,2_shiftH.png + :scale: 25% +.. |f21| image:: images/kinetics_families/Cyclic_Ether_Formation.png + :scale: 25% +.. |f22| image:: images/kinetics_families/Cyclic_Thioether_Formation.png + :scale: 25% +.. |f23| image:: images/kinetics_families/Cyclopentadiene_scission.png + :scale: 25% +.. |f24| image:: images/kinetics_families/Diels_alder_addition.png + :scale: 25% +.. |f25| image:: images/kinetics_families/Disproportionation.png + :scale: 25% +.. |f26| image:: images/kinetics_families/H_Abstraction.png + :scale: 25% +.. |f27| image:: images/kinetics_families/HO2_Elimination_from_PeroxyRadical.png + :scale: 25% +.. |f28| image:: images/kinetics_families/Intra_2+2_cycloaddition_Cd.png + :scale: 25% +.. |f29| image:: images/kinetics_families/Intra_5_membered_conjugated_C=C_C=C_addition.png + :scale: 25% +.. |f30| image:: images/kinetics_families/Intra_Diels_alder_monocyclic.png + :scale: 25% +.. |f31| image:: images/kinetics_families/Intra_Disproportionation.png + :scale: 25% +.. |f32| image:: images/kinetics_families/Intra_ene_reaction.png + :scale: 25% +.. |f33| image:: images/kinetics_families/intra_H_migration.png + :scale: 25% +.. |f34| image:: images/kinetics_families/intra_NO2_ONO_conversion.png + :scale: 25% +.. |f35| image:: images/kinetics_families/intra_OH_migration.png + :scale: 25% +.. |f36| image:: images/kinetics_families/Intra_R_Add_Endocyclic.png + :scale: 25% +.. |f37| image:: images/kinetics_families/Intra_R_Add_Exocyclic.png + :scale: 25% +.. |f39| image:: images/kinetics_families/Intra_R_Add_ExoTetCyclic.png + :scale: 25% +.. |f40| image:: images/kinetics_families/Intra_Retro_Diels_alder_bicyclic.png + :scale: 25% +.. |f41| image:: images/kinetics_families/Intra_RH_Add_Endocyclic.png + :scale: 25% +.. |f42| image:: images/kinetics_families/Intra_RH_Add_Exocyclic.png + :scale: 25% +.. |f43| image:: images/kinetics_families/intra_substitutionCS_cyclization.png + :scale: 25% +.. |f44| image:: images/kinetics_families/intra_substitutionCS_isomerization.png + :scale: 25% +.. |f45| image:: images/kinetics_families/intra_substitutionS_cyclization.png + :scale: 25% +.. |f46| image:: images/kinetics_families/intra_substitutionS_isomerization.png + :scale: 25% +.. |f47| image:: images/kinetics_families/ketoenol.png + :scale: 25% +.. |f48| image:: images/kinetics_families/Korcek_step1.png + :scale: 25% +.. |f49| image:: images/kinetics_families/Korcek_step2.png + :scale: 25% +.. |f50| image:: images/kinetics_families/lone_electron_pair_bond.png + :scale: 25% +.. |f51| image:: images/kinetics_families/R_Addition_COm.png + :scale: 25% +.. |f52| image:: images/kinetics_families/R_Addition_CSm.png + :scale: 25% +.. |f53| image:: images/kinetics_families/R_Addition_MultipleBond.png + :scale: 25% +.. |f54| image:: images/kinetics_families/R_Recombination.png + :scale: 25% +.. |f55| image:: images/kinetics_families/Singlet_Carbene_Intra_Disproportionation.png + :scale: 25% +.. |f56| image:: images/kinetics_families/Singlet_Val6_to_triplet.png + :scale: 25% +.. |f57| image:: images/kinetics_families/Substitution_O.png + :scale: 25% +.. |f58| image:: images/kinetics_families/SubstitutionS.png + :scale: 25% + +.. table:: + :widths: 50 50 + + ===================================================== ===== + **1,2-Birad_to_alkene** |f01| + **1,2_Insertion_carbene** |f02| + **1,2_Insertion_CO** |f03| + **1,2_shiftS** |f05| + **1,3_Insertion_CO2** |f06| + **1,3_Insertion_ROR** |f07| + **1,3_Insertion_RSR** |f08| + **1,4_Cyclic_birad_scission** |f09| + **1,4_Linear_birad_scission** |f10| + **1+2_Cycloaddition** |f11| + **2+2_cycloaddition_CCO** |f12| + **2+2_cycloaddition_Cd** |f13| + **2+2_cycloaddition_CO** |f14| + **2+2_cycloaddition_CS** |f15| + **6_membered_central_C-C_shift** |f16| + **Birad_recombination** |f17| + **Birad_R_Recombination** |f18| + **CO_Disproportionation** |f19| + **Concerted_Intra_Diels_alder_monocyclic_1,2_shiftH** |f20| + **Cyclic_Ether_Formation** |f21| + **Cyclic_Thioether_Formation** |f22| + **Cyclopentadiene_scission** |f23| + **Diels_alder_addition** |f24| + **Disproportionation** |f25| + **H_Abstraction** |f26| + **HO2_Elimination_from_PeroxyRadical** |f27| + **Intra_2+2_cycloaddition_Cd** |f28| + **Intra_5_membered_conjugated_C=C_C=C_addition** |f29| + **Intra_Diels_alder_monocyclic** |f30| + **Intra_Disproportionation** |f31| + **Intra_ene_reaction** |f32| + **intra_H_migration** |f33| + **intra_NO2_ONO_conversion** |f34| + **intra_OH_migration** |f35| + **Intra_R_Add_Endocyclic** |f36| + **Intra_R_Add_Exocyclic** |f37| + **Intra_R_Add_ExoTetCyclic** |f39| + **Intra_Retro_Diels_alder_bicyclic** |f40| + **Intra_RH_Add_Endocyclic** |f41| + **Intra_RH_Add_Exocyclic** |f42| + **intra_substitutionCS_cyclization** |f43| + **intra_substitutionCS_isomerization** |f44| + **intra_substitutionS_cyclization** |f45| + **intra_substitutionS_isomerization** |f46| + **ketoenol** |f47| + **Korcek_step1** |f48| + **Korcek_step2** |f49| + **lone_electron_pair_bond** |f50| + **R_Addition_COm** |f51| + **R_Addition_CSm** |f52| + **R_Addition_MultipleBond** |f53| + **R_Recombination** |f54| + **Singlet_Carbene_Intra_Disproportionation** |f55| + **Singlet_Val6_to_triplet** |f56| + **Substitution_O** |f57| + **SubstitutionS** |f58| + ===================================================== ===== Recipe diff --git a/documentation/source/users/rmg/input.rst b/documentation/source/users/rmg/input.rst index 50de5173a4..d70e504f19 100644 --- a/documentation/source/users/rmg/input.rst +++ b/documentation/source/users/rmg/input.rst @@ -7,15 +7,15 @@ Creating Input Files The syntax and parameters within an RMG input file are explained below. We recommend trying to build your first input file while referencing one of the -:ref:`Example Input Files`. Alternatively, you can use our web form found +:ref:`Example Input Files`. Alternatively, you can use our web form found at http://rmg.mit.edu/input to assist in creating an input file. Syntax ====== -The format of RMG-Py :file:`input.py` is based on Python syntax. +The format of RMG-Py :file:`input.py` is based on Python syntax. -Each section is made up of one or more function calls, where parameters are +Each section is made up of one or more function calls, where parameters are specified as text strings, numbers, or objects. Text strings must be wrapped in either single or double quotes. @@ -35,7 +35,7 @@ you may specify the thermodynamic properties of species in the ThermoLibrary. When a species is specified in the ThermoLibrary, RMG will automatically use those thermodynamic properties instead of generating them from Benson's formulas. Multiple libraries may be created, if so desired. -The order in which the thermo libraries are specified is important: +The order in which the thermo libraries are specified is important: If a species appears in multiple thermo libraries, the first instance will be used. @@ -51,16 +51,16 @@ by Benson's method. For example, if you wish to use the GRI-Mech 3.0 mechanism [GRIMech3.0]_ as a ThermoLibrary in your model, the syntax will be:: thermoLibraries = ['primaryThermoLibrary','GRI-Mech3.0'] - + .. [GRIMech3.0] Gregory P. Smith, David M. Golden, Michael Frenklach, Nigel W. Moriarty, Boris Eiteneer, Mikhail Goldenberg, C. Thomas Bowman, Ronald K. Hanson, Soonho Song, William C. Gardiner, Jr., Vitali V. Lissianski, and Zhiwei Qin http://www.me.berkeley.edu/gri_mech/ -This library is located in the +This library is located in the :file:`$RMG/RMG-database/input/thermo/libraries` directory. All "Locations" for the ThermoLibrary field must be with respect to the :file:`$RMG/RMG-database/input/thermo/libraries` directory. .. note:: - Checks during the initialization are maid to avoid users to use "liquid thermo librairies" in gas phase simulations or to use + Checks during the initialization are maid to avoid users to use "liquid thermo librairies" in gas phase simulations or to use "liquid phase libraries" obtained in another solvent that the one defined in the input file in liquid phase simulations. .. _reactionlibraries: @@ -69,18 +69,18 @@ Reaction Libraries ------------------ The next section of the :file:`input.py` file specifies which, if any, Reaction Libraries should be used. When a reaction library is specified, RMG will first -use the reaction library to generate all the relevant reactions for the species -in the core before going through the reaction templates. Unlike the Seed Mechanism, -reactions present in a Reaction Library will not be included in the core automatically -from the start. - -You can specify your own reaction library in the location section. -In the following example, the user has created -a reaction library with a few additional reactions specific to n-butane, and these reactions +use the reaction library to generate all the relevant reactions for the species +in the core before going through the reaction templates. Unlike the Seed Mechanism, +reactions present in a Reaction Library will not be included in the core automatically +from the start. + +You can specify your own reaction library in the location section. +In the following example, the user has created +a reaction library with a few additional reactions specific to n-butane, and these reactions are to be used in addition to the Glarborg C3 library:: reactionLibraries = [('Glarborg/C3',False)], - + The keyword False/True permits user to append all unused reactions (= kept in the edge) from this library to the chemkin file. True means those reactions will be appended. Using just the string inputs would lead to a default value of `False`. In the previous example, this would look like:: @@ -94,8 +94,8 @@ Because the units for the Arrhenius parameters are given in each mechanism, the different mechanisms can have different units. .. note:: - While using a Reaction Library the user must be careful enough to provide - all instances of a particular reaction in the library file, as RMG will + While using a Reaction Library the user must be careful enough to provide + all instances of a particular reaction in the library file, as RMG will ignore all reactions generated by its templates. For example, suppose you supply the Reaction Library with butyl_1 --> butyl_2. Although RMG would find two unique instances of this reaction (via a three- and four-member cyclic Transition State), @@ -103,13 +103,13 @@ given in each mechanism, the different mechanisms can have different units. RMG will not handle irreversible reactions correctly, if supplied in a Reaction Library. - + .. _seedmechanism: Seed Mechanisms --------------- -The next section of the :file:`input.py` file specifies which, if any, +The next section of the :file:`input.py` file specifies which, if any, Seed Mechanisms should be used. If a seed mechanism is passed to RMG, every species and reaction present in the seed mechanism will be placed into the core, in addition to the species that are listed in the :ref:`species` section. @@ -122,14 +122,14 @@ seed mechanism in the location section. Please note that the oxidation library should not be used for pyrolysis models. The syntax for the seed mechanisms is similar to that of the primary reaction libraries. :: - seedMechanisms = ['GRI-Mech3.0'] + seedMechanisms = ['GRI-Mech3.0'] The seed mechanisms are stored in :file:`RMG-database/input/kinetics/libraries/` -As the units for the Arrhenius parameters are given in each mechanism, -different mechanisms can have different units. Additionally, if the same -reaction occurs more than once in the combined mechanism, +As the units for the Arrhenius parameters are given in each mechanism, +different mechanisms can have different units. Additionally, if the same +reaction occurs more than once in the combined mechanism, the instance of it from the first mechanism in which it appears is the one that gets used. @@ -137,28 +137,28 @@ the one that gets used. Kinetics Depositories --------------------- -:: +:: kineticsDepositories = ['training'] - - + + .. _kineticsfamilies: Kinetics Families ----------------- -In this section users can specify the particular reaction families that they wish to use to generate their model. for example you can use only :file:`Intra_RH_Add_Endocyclic` family to build the model by:: +In this section users can specify the particular reaction families that they wish to use to generate their model. for example you can use only :file:`Intra_RH_Add_Endocyclic` family to build the model by:: kineticsFamilies = ['Intra_RH_Add_Endocyclic'] - + Otherwise, by typing 'default' (and excluding the brackets that are shown in the example above), RMG will use recommended reaction families to generate the mechanism. The recommended reaction families can be found in :file:`RMG-database/input/families/recommended.py`. - + Kinetics Estimator ------------------ -The last section is specifying that RMG is estimating kinetics of reactions from rate rules. For more details on how kinetic estimations is working check :ref:`Kinetics Estimation `:: +The last section is specifying that RMG is estimating kinetics of reactions from rate rules. For more details on how kinetic estimations is working check :ref:`Kinetics Estimation `:: kineticsEstimator = 'rate rules' - + The following is an example of a database block, based on above chosen libraries and options:: @@ -166,7 +166,7 @@ The following is an example of a database block, based on above chosen libraries thermoLibraries = ['primaryThermoLibrary', 'GRI-Mech3.0'], reactionLibraries = [('Glarborg/C3',False)], seedMechanisms = ['GRI-Mech3.0'], - kineticsDepositories = ['training'], + kineticsDepositories = ['training'], kineticsFamilies = 'defult', kineticsEstimator = 'rate rules', ) @@ -176,15 +176,15 @@ The following is an example of a database block, based on above chosen libraries List of species =============== -Species to be included in the core at the start of your RMG job are defined in the species block. +Species to be included in the core at the start of your RMG job are defined in the species block. The label, reactive or inert, and structure of each reactant must be specified. -The label field will be used throughout your mechanism to identify the species. -Inert species in the model can be defined by setting reactive to be ``False``. Reaction -families will no longer be applied to these species, but reactions of the inert from libraries -and seed mechanisms will still be considered. For all other species the reactive status must -be set as ``True``. The structure of the species can be defined using either by using SMILES or -:ref:`adjacencyList `. +The label field will be used throughout your mechanism to identify the species. +Inert species in the model can be defined by setting reactive to be ``False``. Reaction +families will no longer be applied to these species, but reactions of the inert from libraries +and seed mechanisms will still be considered. For all other species the reactive status must +be set as ``True``. The structure of the species can be defined using either by using SMILES or +:ref:`adjacencyList `. The following is an example of a typical species item, based on methane using SMILE or adjacency list to define the structure:: @@ -193,7 +193,7 @@ The following is an example of a typical species item, based on methane using SM reactive=True, structure=SMILES("C"), ) - + species( label='CH4', reactive=True, @@ -210,8 +210,8 @@ Reaction System =============== Every reaction system we want the model to be generated at must be defined individually. -Currently, RMG can only model constant temperature and pressure systems. Future versions -will allow for variable temperature and pressure. To define a reaction system we need to +Currently, RMG can only model constant temperature and pressure systems. Future versions +will allow for variable temperature and pressure. To define a reaction system we need to define the temperature, pressure and initial mole fractions of the reactant species. The initial mole fractions are defined using the label for the species in the species block. Every reaction system can have its termination criterion based on @@ -236,23 +236,23 @@ The following is an example of a simple reactor system:: sensitivityThreshold=0.001, ) - + Troubleshooting tip: if you are using a goal conversion rather than time, the reaction systems may reach equilibrium below the goal conversion, leading to a job that cannot converge physically. Therefore it is may be necessary to reduce the goal conversion or set a goal reaction time. -For sensitivity analysis, RMG-Py must be compiled with the DASPK solver, which is done by default but has -some dependency restrictions. (See :ref:`License Restrictions on Dependencies ` for more details.) +For sensitivity analysis, RMG-Py must be compiled with the DASPK solver, which is done by default but has +some dependency restrictions. (See :ref:`License Restrictions on Dependencies ` for more details.) The sensitivity and sensitivityThrehold are optional arguments for when the user would like to conduct sensitivity analysis with respect to the reaction rate -coefficients for the list of species given for ``sensitivity``. +coefficients for the list of species given for ``sensitivity``. -Sensitivity analysis is conducted for the list of species given for ``sensitivity`` argument in the input file. -The normalized concentration sensitivities with respect to the reaction rate coefficients dln(C_i)/dln(k_j) are saved to a csv file -with the file name ``sensitivity_1_SPC_1.csv`` with the first index value indicating the reactor system and the second naming the index of the species +Sensitivity analysis is conducted for the list of species given for ``sensitivity`` argument in the input file. +The normalized concentration sensitivities with respect to the reaction rate coefficients dln(C_i)/dln(k_j) are saved to a csv file +with the file name ``sensitivity_1_SPC_1.csv`` with the first index value indicating the reactor system and the second naming the index of the species the sensitivity analysis is conducted for. Sensitivities to thermo of individual species is also saved as semi normalized sensitivities dln(C_i)/d(G_j) where the units are given in 1/(kcal mol-1). The sensitivityThreshold is set to some value so that only -sensitivities for dln(C_i)/dln(k_j) > sensitivityThreshold or dlnC_i/d(G_j) > sensitivityThreshold are saved to this file. +sensitivities for dln(C_i)/dln(k_j) > sensitivityThreshold or dlnC_i/d(G_j) > sensitivityThreshold are saved to this file. Note that in the RMG job, after the model has been generated to completion, sensitivity analysis will be conducted in one final simulation (sensitivity is not performed in intermediate iterations of the job). @@ -262,7 +262,7 @@ in one final simulation (sensitivity is not performed in intermediate iterations Simulator Tolerances ==================== The next two lines specify the absolute and relative tolerance for the ODE solver, respectively. Common values for the absolute tolerance are 1e-15 to 1e-25. Relative tolerance is usually 1e-4 to 1e-8:: - + simulator( atol=1e-16, rtol=1e-8, @@ -279,25 +279,25 @@ are set to a default value of 1e-6 and 1e-4 respectively unless the user specifi Model Tolerances ================ -Model tolerances dictate how species get included in the model. For more information, see the theory behind how RMG builds models using the :ref:`Flux-based Algorithm `. +Model tolerances dictate how species get included in the model. For more information, see the theory behind how RMG builds models using the :ref:`Flux-based Algorithm `. For running an initial job, it is recommended to only change the ``toleranceMoveToCore`` and ``toleranceInterruptSimulation`` values to an equivalent desired value. We find that typically a value between ``0.01`` and ``0.05`` is best. If your model cannot converge within a few hours, more advanced settings such as :ref:`reaction filtering ` or :ref:`pruning ` can be turned on to speed up your simulation at a slight risk of omitting chemistry. :: - + model( toleranceMoveToCore=0.1, toleranceInterruptSimulation=0.1, ) -- ``toleranceMoveToCore`` indicates how high the edge flux ratio for a species must get to enter the core model. This tolerance is designed for controlling the accuracy of final model. +- ``toleranceMoveToCore`` indicates how high the edge flux ratio for a species must get to enter the core model. This tolerance is designed for controlling the accuracy of final model. - ``toleranceInterruptSimulation`` indicates how high the edge flux ratio must get to interrupt the simulation (before reaching the ``terminationConversion`` or ``terminationTime``). This value should be set to be equal to ``toleranceMoveToCore`` unless the advanced :ref:`pruning ` feature is desired. .. _filterReactions: Advanced Setting: Speed Up by Filtering Reactions ------------------------------------------------- -For generating models for larger molecules, RMG-Py may have trouble converging because it must find reactions on the order of -:math:`(n_{reaction\: sites})^{{n_{species}}}`. Thus it can be further sped up by pre-filtering reactions that are +For generating models for larger molecules, RMG-Py may have trouble converging because it must find reactions on the order of +:math:`(n_{reaction\: sites})^{{n_{species}}}`. Thus it can be further sped up by pre-filtering reactions that are added to the model. This modification to the algorithm does not react core species together until their concentrations are deemed high enough. It is recommended to turn on this flag when the model does not converge with normal parameter settings. See :ref:`Filtering Reactions within the Flux-based Algorithm `. for more details. :: @@ -307,7 +307,7 @@ the model does not converge with normal parameter settings. See :ref:`Filtering toleranceInterruptSimulation=0.1, filterReactions=True, ) - + **Additional parameters:** - ``filterReactions``: set to ``True`` if reaction filtering is turned on. By default it is set to False. @@ -316,7 +316,7 @@ the model does not converge with normal parameter settings. See :ref:`Filtering Advanced Setting: Speed Up by Pruning ------------------------------------- -For further speed-up, it is also possible to perform mechanism generation with pruning of “unimportant” edge species to reduce memory usage. +For further speed-up, it is also possible to perform mechanism generation with pruning of “unimportant” edge species to reduce memory usage. A typical set of parameters for pruning is:: @@ -324,7 +324,7 @@ A typical set of parameters for pruning is:: toleranceMoveToCore=0.5, toleranceInterruptSimulation=1e8, toleranceKeepInEdge=0.05, - maximumEdgeSpecies=200000 + maximumEdgeSpecies=200000, minCoreSizeForPrune=50, minSpeciesExistIterationsForPrune=2, ) @@ -332,24 +332,24 @@ A typical set of parameters for pruning is:: **Additional parameters:** - ``toleranceKeepInEdge`` indicates how low the edge flux ratio for a species must be to keep on the edge. This should be set to zero, which is its default. -- ``maximumEdgeSpecies`` indicates the upper limit for the size of the edge. The default value is set to ``1000000`` species. -- ``minCoreSizeForPrune`` ensures that a minimum number of species are in the core before pruning occurs, in order to avoid pruning the model when it is far away from completeness. The default value is set to 50 species. -- ``minSpeciesExistIterationsForPrune`` is set so that the edge species stays in the job for at least that many iterations before it can be pruned. The default value is 2 iterations. +- ``maximumEdgeSpecies`` indicates the upper limit for the size of the edge. The default value is set to ``1000000`` species. +- ``minCoreSizeForPrune`` ensures that a minimum number of species are in the core before pruning occurs, in order to avoid pruning the model when it is far away from completeness. The default value is set to 50 species. +- ``minSpeciesExistIterationsForPrune`` is set so that the edge species stays in the job for at least that many iterations before it can be pruned. The default value is 2 iterations. **Recommendations:** We recommend setting ``toleranceKeepInEdge`` to not be larger than 10% of ``toleranceMoveToCore``, based on a pruning case study. -In order to always enable pruning, ``toleranceInterruptSimulation`` should be set as a high value, e.g. 1e8. +In order to always enable pruning, ``toleranceInterruptSimulation`` should be set as a high value, e.g. 1e8. ``maximumEdgeSpecies`` can be adjusted based on user's RAM size. Usually 200000 edge species would cause memory shortage of 8GB computer, setting ``maximumEdgeSpecies = 200000`` (or lower values) could effectively prevent memory crash. **Additional Notes:** -Note that when using pruning, RMG will not prune unless all reaction systems reach the goal reaction time or conversion without exceeding the ``toleranceInterruptSimulation``. +Note that when using pruning, RMG will not prune unless all reaction systems reach the goal reaction time or conversion without exceeding the ``toleranceInterruptSimulation``. Therefore, you may find that RMG is not pruning even though the model edge size exceeds ``maximumEdgeSpecies``, or an edge species has flux below the ``toleranceKeepInEdge``. This is -a safety check within RMG to ensure that species are not pruned too early, resulting in inaccurate chemistry. In order to increase the likelihood of pruning you can -try increasing ``toleranceInterruptSimulation`` to an arbitrarily high value. +a safety check within RMG to ensure that species are not pruned too early, resulting in inaccurate chemistry. In order to increase the likelihood of pruning you can +try increasing ``toleranceInterruptSimulation`` to an arbitrarily high value. As a contrast, a typical set of parameters for non-pruning is:: @@ -359,36 +359,133 @@ As a contrast, a typical set of parameters for non-pruning is:: toleranceInterruptSimulation=0.5, ) -where ``toleranceKeepInEdge`` is always 0, meaning all the edge species will be kept in edge since all the edge species have positive flux. -``toleranceInterruptSimulation`` equals to ``toleranceMoveToCore`` so that ODE simulation get interrupted once discovering a new core species. +where ``toleranceKeepInEdge`` is always 0, meaning all the edge species will be kept in edge since all the edge species have positive flux. +``toleranceInterruptSimulation`` equals to ``toleranceMoveToCore`` so that ODE simulation get interrupted once discovering a new core species. Because the ODE simulation is always interrupted, no pruning is performed. Please find more details about the theory behind pruning at :ref:`Pruning Theory `. +Advanced Setting: Thermodynamic Pruning +---------------------------------------------------- +Thermodynamic pruning is an alternative to flux pruning that does not require a given +simulation to complete to remove excess species. The thermodynamic criteria is +calculated by determining the minimum and maximum Gibbs energies of formation (Gmin and Gmax) +among species in the core. If the Gibbs energy of formation of a given species is G +the value of the criteria is (G-Gmax)/(Gmax-Gmin). All of the Gibbs energies are evaluated +at the highest temperature used in all of the reactor systems. This means that a value of 0.2 +for the criterion implies that it will not add species that have Gibbs energies of formation +greater than 20% of the core Gibbs energy range greater than the maximum Gibbs energy of +formation within the core. + +For example :: + + model( + toleranceMoveToCore=0.5, + toleranceInterruptSimulation=0.5, + toleranceThermoKeepSpeciesInEdge=0.5, + maximumEdgeSpecies=200000, + minCoreSizeForPrune=50, + ) + +**Advantages over flux pruning**: + +Species are removed immediately if they violate tolerance +Completing a simulation is unnecessary for this pruning so there is no need +to waste time setting the interrupt tolerance higher than the movement tolerance. +Will always maintain the correct ``maximumEdgeSpecies``. + +**Primary disadvantage**: + +Since we determine whether to add species primarily based on flux, at tight tolerances +this is more likely to kick out species RMG might otherwise have added to core. + + Advanced Setting: Taking Multiple Species At A Time ---------------------------------------------------- Taking multiple objects (species, reactions or pdepNetworks) during a given simulation can often decrease your overall model generation time -over only taking one. For this purpose there is a maxNumObjsPerIter parameter that allows RMG to take -that many species, reactions or pdepNetworks from a given simulation. This is done in the order they trigger their respective criteria. +over only taking one. For this purpose there is a ``maxNumObjsPerIter`` parameter that allows RMG to take +that many species, reactions or pdepNetworks from a given simulation. This is done in the order they trigger their respective criteria. + +You can also set ``terminateAtMaxObjects=True`` to cause it to terminate when it has the maximum +number of objects allowed rather than waiting around until it hits an interrupt tolerance. This +avoids additional simulation time, but will also make it less likely to finish simulations, which can +affect flux pruning. For example :: model( toleranceKeepInEdge=0.0, toleranceMoveToCore=0.1, - toleranceInterruptSimulation=0.1, + toleranceInterruptSimulation=0.3, maxNumObjsPerIter=2, + terminateAtMaxObjects=True, ) Note that this can also result in larger models, however, sometimes these larger models (from taking more than one -object at a time) pick up chemistry that would otherwise have been missed. +object at a time) pick up chemistry that would otherwise have been missed. .. _ontheflyquantumcalculations: +Advanced Setting: Dynamics Criterion +------------------------------------- +While the flux criterion works very well for identifying new species that have high flux +and therefore will likely be high throughput or high concentration species, +it provides few automatic guarantees about how well a given model will accurately represent +the concentrations of the involved species. The dynamics criterion is more complex +than the flux criterion, but in general it is a measure of how much impact a given +reaction will have on the concentrations of core species. A more detailed explanation +is available in the theory section: :ref:`dynamics`. + +Reasonable values for the dynamics criterion range typically between 2-30. + +For example :: + + model( + toleranceMoveToCore=0.1, + toleranceInterruptSimulation=0.1, + toleranceMoveEdgeReactionToCore=10.0, + toleranceMoveEdgeReactionToCoreInterrupt=5.0, + ) + +Note that it is highly recommended to use the dynamics criterion only alongside the flux +criterion and not by itself. + +Advanced Setting: Surface Algorithm +------------------------------------- +One common issue with the dynamics criterion is that it treats all core species equally +as discussed in our theory section: :ref:`dynamics`. Because of this, if the dynamics +criterion is set too low it enters a feedback loop where it adds species and then +since it can't get those species' concentrations right it adds more species and so on. +In order to avoid this feedback loop the surface algorithm was developed. It creates +a new partition called the *surface* that is considered part of the core. We will +refer to the part of the core that is not part of the surface as the *bulk core*. When +operating without the dynamics criterion everything moves from edge to the bulk core as usual; +however the dynamics criterion is managed differently. When using the surface algorithm most +reactions pulled in by the dynamics criterion enter the surface instead of the bulk core. +However, unlike movement to bulk core a constraint is placed on movement to the surface. +Any reaction moved to the surface must have either both reactants or both products +in the bulk core. This prevents the dynamics criterion from pulling in reactions +to get the concentrations of species in the surface right avoiding the feedback loop. +To avoid important species being trapped in the surface we also add criteria +for movement from surface to bulk core based on flux or dynamics criterion. + + +For example :: + + model( + toleranceMoveToCore=0.1, + toleranceInterruptSimulation=0.1, + toleranceMoveEdgeReactionToCore=30.0, + toleranceMoveEdgeReactionToCoreInterrupt=5.0, + toleranceMoveEdgeReactionToSurface=10.0, + toleranceMoveSurfaceSpeciesToCore=.01, + toleranceMoveSurfaceReactionToCore=5.0, + ) + On the fly Quantum Calculations =============================== -This block is used when quantum mechanical calculations are desired to determine thermodynamic parameters. +This block is used when quantum mechanical calculations are desired to determine thermodynamic parameters. These calculations are only run if the molecule is not included in a specified thermo library. The ``onlyCyclics`` option, if ``True``, only runs these calculations for cyclic species. In this case, group additive estimates are used for all other species. @@ -401,7 +498,7 @@ A folder can be specified to store the files used in these calculations, however if not specified this defaults to a `QMfiles` folder in the output folder. The calculations are also only run on species with a maximum radical number set by the user. -If a molecule has a higher radical number, the molecule is saturated with hydrogen atoms, then +If a molecule has a higher radical number, the molecule is saturated with hydrogen atoms, then quantum mechanical calculations with subsequent hydrogen bond incrementation is used to determine the thermodynamic parameters. @@ -424,11 +521,11 @@ The following is an example of the quantum mechanics options :: Pressure Dependence =================== -This block is used when the model should account for pressure -dependent rate coefficients. RMG can estimate pressure dependence kinetics based on ``Modified Strong Collision`` and ``Reservoir State`` methods. -The former utilizes the modified strong collision approach of Chang, Bozzelli, and Dean [Chang2000]_, -and works reasonably well while running more rapidly. The latter -utilizes the steady-state/reservoir-state approach of Green and Bhatti [Green2007]_, +This block is used when the model should account for pressure +dependent rate coefficients. RMG can estimate pressure dependence kinetics based on ``Modified Strong Collision`` and ``Reservoir State`` methods. +The former utilizes the modified strong collision approach of Chang, Bozzelli, and Dean [Chang2000]_, +and works reasonably well while running more rapidly. The latter +utilizes the steady-state/reservoir-state approach of Green and Bhatti [Green2007]_, and is more theoretically sound but more expensive. @@ -479,11 +576,11 @@ Temperature and pressure for the interpolation scheme To generate the :math:`k(T,P)` interpolation model, a set of temperatures and pressures must be used. RMG can do this automatically, but it must be told a few parameters. We need to specify the limits of the temperature and pressure for the fitting of the interpolation scheme and the number of points to be considered in between this limit. -For typical combustion model temperatures of the experiments range from 300 - 2000 K and pressure 1E-2 to 100 bar :: +For typical combustion model temperatures of the experiments range from 300 - 2000 K and pressure 1E-2 to 100 bar :: temperatures=(300,2000,'K',8) pressures=(0.01,100,'bar',5) - + Interpolation scheme -------------------- @@ -491,27 +588,27 @@ To use logarithmic interpolation of pressure and Arrhenius interpolation for tem line :: interpolation=('PDepArrhenius',) - + The auxillary information printed to the Chemkin chem.inp file will have the "PLOG" -format. Refer to Section 3.5.3 of the :file:`CHEMKIN_Input.pdf` document and/or +format. Refer to Section 3.5.3 of the :file:`CHEMKIN_Input.pdf` document and/or Section 3.6.3 of the :file:`CHEMKIN_Theory.pdf` document. These files are part of -the CHEMKIN manual. +the CHEMKIN manual. -To fit a set of Chebyshev polynomials on inverse temperature and logarithmic pressure axes mapped -to [-1,1], specify `''Chebyshev'` interpolation. -You should also specify the number of temperature and pressure basis functions by adding the appropriate integers. +To fit a set of Chebyshev polynomials on inverse temperature and logarithmic pressure axes mapped +to [-1,1], specify `''Chebyshev'` interpolation. +You should also specify the number of temperature and pressure basis functions by adding the appropriate integers. For example, the following specifies that six basis functions in temperature and four in pressure should be used :: interpolation=('Chebyshev', 6, 4) The auxillary information printed to the Chemkin chem.inp file will have the "CHEB" -format. Refer to Section 3.5.3 of the :file:`CHEMKIN_Input.pdf` document and/or +format. Refer to Section 3.5.3 of the :file:`CHEMKIN_Input.pdf` document and/or Section 3.6.4 of the :file:`CHEMKIN_Theory.pdf` document. Regarding the number of polynomial coeffients for Chebyshev interpolated rates, -plese refer to the :class:`rmgpy.kinetics.Chebyshev` documentation. -The number of pressures and temperature coefficents should always be smaller -than the respective number of user-specified temperatures and pressures. +plese refer to the :class:`rmgpy.kinetics.Chebyshev` documentation. +The number of pressures and temperature coefficents should always be smaller +than the respective number of user-specified temperatures and pressures. Maximum size of adduct for which pressure dependence kinetics be generated @@ -520,7 +617,7 @@ Maximum size of adduct for which pressure dependence kinetics be generated By default pressure dependence is run for every system that might show pressure dependence, i.e. every isomerization, dissociation, and association reaction. In reality, larger molecules are less likely to exhibit pressure-dependent -behavior than smaller molecules due to the presence of more modes for +behavior than smaller molecules due to the presence of more modes for randomization of the internal energy. In certain cases involving very large molecules, it makes sense to only consider pressure dependence for molecules smaller than some user-defined number of atoms. This is specified e.g. using @@ -535,9 +632,9 @@ of atoms (16 in the above example). .. _miscellaneousoptions: Miscellaneous Options -===================== +===================== -Miscellaneous options:: +Miscellaneous options:: options( name='Seed', @@ -563,21 +660,21 @@ The ``units`` field is set to ``si``. Currently there are no other unit options The ``saveRestartPeriod`` indictes how frequently you wish to save restart files. For very large/long RMG jobs, this process can take a significant amount of time. In such cases, the user may wish to increase the time period for saving these restart files. Setting ``generateOutputHTML`` to ``True`` will let RMG know that you want to save 2-D images (png files in the local ``species`` folder) of all species in the generated core model. It will save a visualized -HTML file for your model containing all the species and reactions. Turning this feature off by setting it to ``False`` may save memory if running large jobs. +HTML file for your model containing all the species and reactions. Turning this feature off by setting it to ``False`` may save memory if running large jobs. Setting ``generatePlots`` to ``True`` will generate a number of plots describing the statistics of the RMG job, including the reaction model core and edge size and memory use versus execution time. These will be placed in the output directory in the plot/ folder. Setting ``saveSimulationProfiles`` to ``True`` will make RMG save csv files of the simulation in .csv files in the ``solver/`` folder. The filename will be ``simulation_1_26.csv`` where the first number corresponds to the reaciton system, and the second number corresponds to the total number of species at the point of the simulation. Therefore, the highest second number will indicate the latest simulation that RMG has complete while enlarging the core model. The information inside the csv file will provide the time, reactor volume in m^3, as well as mole fractions of the individual species. -Setting ``verboseComments`` to ``True`` will make RMG generate chemkin files with complete verbose commentary for the kinetic and thermo parameters. This will be helpful in debugging what values are being averaged for the kinetics. Note that this may produce very large files. +Setting ``verboseComments`` to ``True`` will make RMG generate chemkin files with complete verbose commentary for the kinetic and thermo parameters. This will be helpful in debugging what values are being averaged for the kinetics. Note that this may produce very large files. -Setting ``saveEdgeSpecies`` to ``True`` will make RMG generate chemkin files of the edge reactions in addition to the core model in files such as ``chem_edge.inp`` and ``chem_edge_annotated.inp`` files located inside the ``chemkin`` folder. These files will be helpful in viewing RMG's estimate for edge reactions and seeing if certain reactions one expects are actually in the edge or not. +Setting ``saveEdgeSpecies`` to ``True`` will make RMG generate chemkin files of the edge reactions in addition to the core model in files such as ``chem_edge.inp`` and ``chem_edge_annotated.inp`` files located inside the ``chemkin`` folder. These files will be helpful in viewing RMG's estimate for edge reactions and seeing if certain reactions one expects are actually in the edge or not. Setting ``keepIrreversible`` to ``True`` will make RMG import library reactions as is, whether they are reversible or irreversible in the library. Otherwise, if ``False`` (default value), RMG will force all library reactions to be reversible, and will assign the forward rate from the relevant library. Species Constraints -===================== +===================== RMG can generate mechanisms with a number of optional species constraints, such as total number of carbon atoms or electrons per species. These are applied to @@ -598,26 +695,26 @@ all of RMG's reaction families. :: allowSingletO2 = False, ) -An additional flag ``allowed`` can be set to allow species +An additional flag ``allowed`` can be set to allow species from either the input file, seed mechanisms, or reaction libraries to bypass these constraints. Note that this should be done with caution, since the constraints will still apply to subsequent -products that form. +products that form. -Note that under all circumstances all forbidden species will still be banned unless they are -manually removed from the database. See :ref:`kineticsDatabase` for more information on -forbidden groups. +Note that under all circumstances all forbidden species will still be banned unless they are +manually removed from the database. See :ref:`kineticsDatabase` for more information on +forbidden groups. -By default, the ``allowSingletO2`` flag is set to ``False``. See :ref:`representing_oxygen` for more information. +By default, the ``allowSingletO2`` flag is set to ``False``. See :ref:`representing_oxygen` for more information. Staging ======== -It is now possible to concatenate different model and simulator blocks into the same run in stages. Any given stage will terminate when the RMG run terminates and then the current group of model and simulator parameters will be switched out with the next group and the run will continue until that stage terminates. Once the last stage terminates the run ends normally. This is currently enabled only for the model and simulator blocks. +It is now possible to concatenate different model and simulator blocks into the same run in stages. Any given stage will terminate when the RMG run terminates and then the current group of model and simulator parameters will be switched out with the next group and the run will continue until that stage terminates. Once the last stage terminates the run ends normally. This is currently enabled only for the model and simulator blocks. -There must be the same number of each of these blocks (although only having one simulator block and many model blocks is enabled as well) and RMG will enter each stage these define in the order they were put in the input file. +There must be the same number of each of these blocks (although only having one simulator block and many model blocks is enabled as well) and RMG will enter each stage these define in the order they were put in the input file. -To enable easier manipulation of staging a new parameter in the model block was developed maxNumSpecies that is the number of core species at which that stage (or if it is the last stage the entire model generation process) will terminate. +To enable easier manipulation of staging a new parameter in the model block was developed maxNumSpecies that is the number of core species at which that stage (or if it is the last stage the entire model generation process) will terminate. For example :: diff --git a/documentation/source/users/rmg/introduction.rst b/documentation/source/users/rmg/introduction.rst index c59a72f743..fcad088b9f 100644 --- a/documentation/source/users/rmg/introduction.rst +++ b/documentation/source/users/rmg/introduction.rst @@ -14,3 +14,4 @@ License RMG is an open source program, available to the general public free of charge. The primary RMG code is distributed under the terms of the `MIT/X11 License `_. However, RMG has a number of dependencies of various licenses, some of which may be more restrictive. **It is the user's responsibility to ensure these licenses have been obtained.** .. literalinclude:: ../../../../LICENSE.txt + :language: none diff --git a/documentation/source/users/rmg/modules/generateFluxDiagram.rst b/documentation/source/users/rmg/modules/generateFluxDiagram.rst index 389dab03a5..ed800302d6 100644 --- a/documentation/source/users/rmg/modules/generateFluxDiagram.rst +++ b/documentation/source/users/rmg/modules/generateFluxDiagram.rst @@ -7,12 +7,32 @@ Generating Flux Diagrams The script, ``generateFluxDiagrams.py``, will create a movie out of a completed RMG model that shows interconnected arrows between species that represent fluxes. -To use this method, you just need a completed RMG run. THe syntax is as follows:: +To use this method, you just need a Chemkin input file and an RMG species dictionary. +The syntax is as follows:: - python generateFluxDiagram.py input.py - -where ``input.py`` is the input file for the completed RMG run. The program will use the automatically -generated file structure to find the other necessary files to create the movie. + python generateFluxDiagram.py [-h] [--java] [--no-dlim] [-s SPECIES] [-f] + [-n N] [-e N] [-c TOL] [-r TOL] [-t S] + INPUT CHEMKIN DICTIONARY [CHEMKIN_OUTPUT] + +Positional arguments:: + + INPUT RMG input file + CHEMKIN Chemkin file + DICTIONARY RMG dictionary file + CHEMKIN_OUTPUT Chemkin output file + +Optional arguments:: + + -h, --help show this help message and exit + --java process RMG-Java model + --no-dlim Turn off diffusion-limited rates + -s DIR, --species DIR Path to folder containing species images + -f, --foreign Not an RMG generated Chemkin file (will be checked for duplicates) + -n N, --maxnode N Maximum number of nodes to show in diagram + -e N, --maxedge N Maximum number of edges to show in diagram + -c TOL, --conctol TOL Lowest fractional concentration to show + -r TOL, --ratetol TOL Lowest fractional species rate to show + -t S, --tstep S Multiplicative factor to use between consecutive time points This method is also available to use with a web browser from the RMG website: `Generate Flux Diagram `_. diff --git a/documentation/source/users/rmg/modules/index.rst b/documentation/source/users/rmg/modules/index.rst index 2655bfd817..e2daabd86a 100644 --- a/documentation/source/users/rmg/modules/index.rst +++ b/documentation/source/users/rmg/modules/index.rst @@ -15,7 +15,7 @@ otherwise. diffModels mergeModels generateReactions - sensitivity + simulate generateFluxDiagram thermoEstimation convertFAME diff --git a/documentation/source/users/rmg/modules/sensitivity.rst b/documentation/source/users/rmg/modules/simulate.rst similarity index 64% rename from documentation/source/users/rmg/modules/sensitivity.rst rename to documentation/source/users/rmg/modules/simulate.rst index 227f6c8db2..ffbeab5397 100644 --- a/documentation/source/users/rmg/modules/sensitivity.rst +++ b/documentation/source/users/rmg/modules/simulate.rst @@ -1,17 +1,17 @@ -.. _sensitivity: +.. _simulate: -******************** -Sensitivity Analysis -******************** +*********************************** +Simulation and Sensitivity Analysis +*********************************** - -For sensitivity analysis, RMG-Py must be compiled with the DASPK solver, which is done by default but has +For sensitivity analysis, RMG-Py must be compiled with the DASPK solver, which is done by default but has some dependency restrictions. (See :ref:`License Restrictions on Dependencies ` for more details.) -Sensitivity analysis can be conducted in a standalone system for an existing kinetics model in Chemkin format. +Sensitivity analysis or a simulation (without sensitivity) can be conducted in a standalone system for an existing +kinetics model in Chemkin format. -To use the sensitivity analysis standalone module:: +To run a simulation and/or sensitivity analysis, use the simulate module in RMG-Py/scripts:: - python sensitivity.py input.py chem.inp species_dictionary.txt + python simulate.py input.py chem.inp species_dictionary.txt where ``chem.inp`` is the CHEMKIN file and the ``species_dictionary.txt`` contains the dictionary of species associated with the CHEMKIN file. ``input.py`` is an input file similar to one used for an RMG job but @@ -22,9 +22,14 @@ does not generate a RMG job. See the following ``input.py`` example file found .. literalinclude:: ../../../../../examples/rmg/minimal_sensitivity/input.py -The names of species named in the input file must coincide with the name specified in the CHEMKIN file. +The names of species named in the input file must coincide with the names specified in the CHEMKIN file. + +Other options that can be specified for the ``simulate.py`` scripts are:: + + --no-dlim Turn off diffusion-limited rates for LiquidReactor + -f, --foreign Not an RMG generated Chemkin file (will be checked for duplicates) -Sensitivity analysis is conducted for the list of species given for ``sensitivity`` argument in the input file. +Sensitivity analysis is conducted for the list of species given for the ``sensitivity`` argument in the input file. The normalized concentration sensitivities with respect to the reaction rate coefficients dln(C_i)/dln(k_j) are saved to a csv file with the file name ``sensitivity_1_SPC_1.csv`` with the first index value indicating the reactor system and the second naming the index of the species the sensitivity analysis is conducted for. Sensitivities to thermo of individual species is also saved as semi normalized sensitivities diff --git a/documentation/source/users/rmg/modules/thermoEstimation.rst b/documentation/source/users/rmg/modules/thermoEstimation.rst index 3a02a43095..c536c4dd99 100644 --- a/documentation/source/users/rmg/modules/thermoEstimation.rst +++ b/documentation/source/users/rmg/modules/thermoEstimation.rst @@ -68,4 +68,8 @@ Submitting a job is easy:: We recommend you make a job-specific directory for each thermoEstimator simulation. +You can also specify that an RMG-style thermo library be saved upon completion:: + + python thermoEstimator.py -l input.py + Note that the RMG website also provides thermo estimation through the `Molecule Search `_. diff --git a/documentation/source/users/rmg/releaseNotes.rst b/documentation/source/users/rmg/releaseNotes.rst index 6873b9da34..9fed704ef2 100644 --- a/documentation/source/users/rmg/releaseNotes.rst +++ b/documentation/source/users/rmg/releaseNotes.rst @@ -4,6 +4,40 @@ Release Notes ************* +RMG-Py Version 2.1.8 +==================== +Date: March 22, 2018 + +- New features: + - Chlorine and iodine atom types have been added, bringing support for these elements to RMG-database + - Forbidden structures now support Molecule and Species definitions in addition to Group definitions + +- Changes: + - Reaction pair generation will now fall back to generic method instead of raising an exception + - Removed sensitivity.py script since it was effectively a duplicate of simulate.py + - Thermo jobs in Cantherm now output a species dictionary + - Fitted atom energy corrections added for B3LYP/6-31g** + - Initial framework added for hydrogen bonding + - Renamed molepro module and associated classes to molpro (MolPro) to match actual spelling of the program + - Chemkin module is now cythonized to improve performance + +- Fixes: + - Allow delocalization of triradicals to prevent hysteresis in resonance structure generation + - Fix reaction comment parsing issue with uncertainty analysis + - Fix numerical issue causing a number of pressure dependent RMG jobs to crash + - Template reactions from seed mechanisms are now loaded as library reactions if the original family is not loaded + - Fix issues with degeneracy calculation for identical reactants + +RMG-database Version 2.1.8 +========================== +Date: March 22, 2018 + +- Changes: + - Corrected name of JetSurf2.0 kinetics and thermo libraries to JetSurf1.0 + - Added actual JetSurf2.0 kinetics and thermo libraries + - Updated thermo groups for near-aromatic radicals, including radical and polycyclic corrections + + RMG-Py Version 2.1.7 ==================== Date: February 12, 2018 diff --git a/examples/cantherm/species/C2H4/ethene.py b/examples/cantherm/species/C2H4/ethene.py index 7fbcbaa41b..f0fd6bb923 100644 --- a/examples/cantherm/species/C2H4/ethene.py +++ b/examples/cantherm/species/C2H4/ethene.py @@ -22,7 +22,7 @@ energy = { 'CBS-QB3': GaussianLog('ethene.log'), 'Klip_2': -78.42735579, - 'CCSD(T)-F12/cc-pVTZ-F12': MoleProLog('ethene_f12.out'), + 'CCSD(T)-F12/cc-pVTZ-F12': MolproLog('ethene_f12.out'), } geometry = GaussianLog('ethene.log') diff --git a/examples/sensitivity/run.sh b/examples/sensitivity/run.sh index 714bf04df3..a0ec3759b6 100755 --- a/examples/sensitivity/run.sh +++ b/examples/sensitivity/run.sh @@ -1,4 +1,4 @@ #!/bin/bash -# Run the stand-alone sensitivity analysis on a kinetic mechanism -python ../../scripts/sensitivity.py input.py chem.inp species_dictionary.txt +# Run the stand-alone sensitivity analysis via simulate.py on a kinetic mechanism +python ../../scripts/simulate.py input.py chem.inp species_dictionary.txt diff --git a/meta.yaml b/meta.yaml index 0e65b76a37..b068c51866 100644 --- a/meta.yaml +++ b/meta.yaml @@ -77,7 +77,7 @@ requirements: - pyzmq - quantities - rdkit >=2015.09.2 - - rmgdatabase >=2.1.7 + - rmgdatabase >=2.1.8 - scipy - scoop - symmetry diff --git a/rmgpy/cantherm/kinetics.py b/rmgpy/cantherm/kinetics.py index 6bd3dae955..e953aa08fa 100644 --- a/rmgpy/cantherm/kinetics.py +++ b/rmgpy/cantherm/kinetics.py @@ -127,6 +127,8 @@ def execute(self, outputFile=None, plot=False): self.save(outputFile) if plot: self.plot(os.path.dirname(outputFile)) + logging.debug('Finished kinetics job for reaction {0}.'.format(self.reaction)) + logging.debug(repr(self.reaction)) def generateKinetics(self,Tlist=None): """ @@ -144,10 +146,12 @@ def generateKinetics(self,Tlist=None): tunneling.E0_TS = (self.reaction.transitionState.conformer.E0.value_si*0.001,"kJ/mol") tunneling.E0_prod = (sum([product.conformer.E0.value_si for product in self.reaction.products])*0.001,"kJ/mol") elif tunneling is not None: - raise ValueError('Unknown tunneling model {0!r}.'.format(tunneling)) - - logging.info('Generating {0} kinetics model for {0}...'.format(kineticsClass, self.reaction)) - + if tunneling.frequency is not None: + # Frequency was given by the user + pass + else: + raise ValueError('Unknown tunneling model {0!r} for reaction {1}.'.format(tunneling, self.reaction)) + logging.debug('Generating {0} kinetics model for {1}...'.format(kineticsClass, self.reaction)) if Tlist is None: Tlist = 1000.0/numpy.arange(0.4, 3.35, 0.05) klist = numpy.zeros_like(Tlist) diff --git a/rmgpy/cantherm/molepro.py b/rmgpy/cantherm/molpro.py similarity index 88% rename from rmgpy/cantherm/molepro.py rename to rmgpy/cantherm/molpro.py index 970a21b801..cd0769dbfd 100644 --- a/rmgpy/cantherm/molepro.py +++ b/rmgpy/cantherm/molpro.py @@ -32,10 +32,10 @@ import rmgpy.constants as constants -class MoleProLog: +class MolproLog: """ - Represents a MolePro log file. The attribute `path` refers to the - location on disk of the MolePro log file of interest. Methods are provided + Represents a Molpro log file. The attribute `path` refers to the + location on disk of the Molpro log file of interest. Methods are provided to extract a variety of information into CanTherm classes and/or NumPy arrays. """ @@ -45,9 +45,9 @@ def __init__(self, path): def loadEnergy(self,frequencyScaleFactor=1.): """ - Return the f12 energy in J/mol from a MolePro Logfile of a CCSD(T)-f12 job. + Return the f12 energy in J/mol from a Molpro Logfile of a CCSD(T)-f12 job. This function determines which energy (f12a or f12b) to use based on the basis set, - which it will parse out of the molepro file. For the vtz and dtz basis sets f12a is + which it will parse out of the Molpro file. For the vtz and dtz basis sets f12a is better approximation, but for higher basis sets f12b is a better approximation """ f = open(self.path, 'r') @@ -61,7 +61,7 @@ def loadEnergy(self,frequencyScaleFactor=1.): else: f12a=False break line=f.readline() - else: raise Exception('Could not find basis set in MolePro File') + else: raise Exception('Could not find basis set in Molpro File') #search for energy E0=None if f12a: @@ -87,4 +87,4 @@ def loadEnergy(self,frequencyScaleFactor=1.): if E0 is not None: E0 = E0 * constants.E_h * constants.Na return E0 - else: raise Exception('Unable to find energy in MolePro log file.') + else: raise Exception('Unable to find energy in Molpro log file.') diff --git a/rmgpy/cantherm/moleproTest.py b/rmgpy/cantherm/molproTest.py similarity index 80% rename from rmgpy/cantherm/moleproTest.py rename to rmgpy/cantherm/molproTest.py index ae382768ac..8aebbc5ad1 100644 --- a/rmgpy/cantherm/moleproTest.py +++ b/rmgpy/cantherm/molproTest.py @@ -32,46 +32,46 @@ import unittest import os -from rmgpy.cantherm.molepro import MoleProLog +from rmgpy.cantherm.molpro import MolproLog import rmgpy.constants as constants ################################################################################ -class MoleProTest(unittest.TestCase): +class MolproTest(unittest.TestCase): """ Contains unit tests for the chempy.io.gaussian module, used for reading - and writing Molepro files. + and writing Molpro files. """ - def testLoadDzFromMoleProLog_F12(self): + def testLoadDzFromMolproLog_F12(self): """ - Uses a Molepro log file for ethylene_dz (C2H4) to test that F12a + Uses a Molpro log file for ethylene_dz (C2H4) to test that F12a energy can be properly read. """ - log=MoleProLog(os.path.join(os.path.dirname(__file__),'data','ethylene_f12_dz.out')) + log=MolproLog(os.path.join(os.path.dirname(__file__),'data','ethylene_f12_dz.out')) E0=log.loadEnergy() self.assertAlmostEqual(E0 / constants.Na / constants.E_h, -78.474353559604, 5) - def testLoadQzFromMoleProLog_F12(self): + def testLoadQzFromMolproLog_F12(self): """ - Uses a Molepro log file for ethylene_qz (C2H4) to test that F12b + Uses a Molpro log file for ethylene_qz (C2H4) to test that F12b energy can be properly read. """ - log=MoleProLog(os.path.join(os.path.dirname(__file__),'data','ethylene_f12_qz.out')) + log=MolproLog(os.path.join(os.path.dirname(__file__),'data','ethylene_f12_qz.out')) E0=log.loadEnergy() self.assertAlmostEqual(E0 / constants.Na / constants.E_h, -78.472682755635, 5) - def testLoadRadFromMoleProLog_F12(self): + def testLoadRadFromMolproLog_F12(self): """ - Uses a Molepro log file for OH (C2H4) to test that radical + Uses a Molpro log file for OH (C2H4) to test that radical energy can be properly read. """ - log=MoleProLog(os.path.join(os.path.dirname(__file__),'data','OH_f12.out')) + log=MolproLog(os.path.join(os.path.dirname(__file__),'data','OH_f12.out')) E0=log.loadEnergy() self.assertAlmostEqual(E0 / constants.Na / constants.E_h, -75.663696424380, 5) diff --git a/rmgpy/cantherm/pdep.py b/rmgpy/cantherm/pdep.py index a84392f112..6c48ebd869 100644 --- a/rmgpy/cantherm/pdep.py +++ b/rmgpy/cantherm/pdep.py @@ -238,6 +238,8 @@ def execute(self, outputFile, plot, format='pdf'): self.save(outputFile) if plot: self.plot(os.path.dirname(outputFile)) + logging.debug('Finished pdep job for reaction {0}.'.format(self.network.label)) + logging.debug(repr(self.network)) def generateTemperatureList(self): """ @@ -362,7 +364,7 @@ def fitInterpolationModels(self): order = len(reaction.reactants) kdata *= 1e6 ** (order-1) kunits = {1: 's^-1', 2: 'cm^3/(mol*s)', 3: 'cm^6/(mol^2*s)'}[order] - + logging.debug('Fitting master eqn data to kinetics for reaction {}.'.format(reaction)) reaction.kinetics = self.fitInterpolationModel(Tdata, Pdata, kdata, kunits) self.network.netReactions.append(reaction) diff --git a/rmgpy/cantherm/statmech.py b/rmgpy/cantherm/statmech.py index 6517a58da3..1881b799a9 100644 --- a/rmgpy/cantherm/statmech.py +++ b/rmgpy/cantherm/statmech.py @@ -43,7 +43,7 @@ from rmgpy.cantherm.output import prettify from rmgpy.cantherm.gaussian import GaussianLog -from rmgpy.cantherm.molepro import MoleProLog +from rmgpy.cantherm.molpro import MolproLog from rmgpy.cantherm.qchem import QchemLog from rmgpy.species import TransitionState, Species @@ -174,6 +174,8 @@ def execute(self, outputFile=None, plot=False): self.load() if outputFile is not None: self.save(outputFile) + logging.debug('Finished statmech job for species {0}.'.format(self.species)) + logging.debug(repr(self.species)) def load(self): """ @@ -200,7 +202,7 @@ def load(self): # File formats 'GaussianLog': GaussianLog, 'QchemLog': QchemLog, - 'MoleProLog': MoleProLog, + 'MolproLog': MolproLog, 'ScanLog': ScanLog, } @@ -262,7 +264,7 @@ def load(self): elif isinstance(energy, QchemLog): energyLog = energy; E0 = None energyLog.path = os.path.join(directory, energyLog.path) - elif isinstance(energy, MoleProLog): + elif isinstance(energy, MolproLog): energyLog = energy; E0 = None energyLog.path = os.path.join(directory, energyLog.path) elif isinstance(energy, float): @@ -646,7 +648,9 @@ def applyEnergyCorrections(E0, modelChemistry, atoms, bonds): elif modelChemistry in ['BMK/cbsb7', 'BMK/6-311G(2d,d,p)']: atomEnergies = {'H':-0.498618853119+ SOC['H'], 'N':-54.5697851544+ SOC['N'], 'O':-75.0515210278+ SOC['O'], 'C':-37.8287310027+ SOC['C'], 'P':-341.167615941+ SOC['P'], 'S': -398.001619915+ SOC['S']} - + elif modelChemistry == 'b3lyp/6-31G**': + atomEnergies = {'H':-0.500426155, 'C':-37.850331697831, 'O':-75.0535872748806, 'S':-398.100820107242} + else: logging.warning('Unknown model chemistry "{0}"; not applying energy corrections.'.format(modelChemistry)) return E0 @@ -701,12 +705,12 @@ def applyEnergyCorrections(E0, modelChemistry, atoms, bonds): 'N-H': 0.06, 'N-N': -0.23, 'N=N': -0.37, 'N#N': -0.64,} elif modelChemistry == 'CBS-QB3': bondEnergies = { - 'C-C': -0.495, 'C-H': -0.045, 'C=C': -0.825, 'C-O': 0.378, 'C=O': 0.743, 'O-H': -0.423, #Table2: Paraskevas, PD (2013). Chemistry-A European J., DOI: 10.1002/chem.201301381 - 'C#C': -0.64, 'C#N': -0.89, 'C-S': 0.43, 'O=S': -0.78, 'S-H': 0.0, 'C-N': -0.13, # Table IX: Petersson GA (1998) J. of Chemical Physics, DOI: 10.1063/1.477794 - 'N-H': -0.42, 'N=O': 1.11, 'N-N': -1.87, 'N=N': -1.58, 'N-O': 0.35, #Table 2: Ashcraft R (2007) J. Phys. Chem. B; DOI: 10.1021/jp073539t - 'N#N': -2.0, 'O=O': -0.2, 'H-H': 1.1, # Unknown source + 'C-C': -0.495,'C-H': -0.045,'C=C': -0.825,'C-O': 0.378,'C=O': 0.743,'O-H': -0.423, #Table2: Paraskevas, PD (2013). Chemistry-A European J., DOI: 10.1002/chem.201301381 + 'C#C': -0.64, 'C#N': -0.89, 'C-S': 0.43, 'O=S': -0.78,'S-H': 0.0, 'C-N': -0.13, 'C-Cl': 1.29, 'C-F': 0.55, # Table IX: Petersson GA (1998) J. of Chemical Physics, DOI: 10.1063/1.477794 + 'N-H': -0.42, 'N=O': 1.11, 'N-N': -1.87, 'N=N': -1.58,'N-O': 0.35, #Table 2: Ashcraft R (2007) J. Phys. Chem. B; DOI: 10.1021/jp073539t + 'N#N': -2.0, 'O=O': -0.2, 'H-H': 1.1, # Unknown source } - elif modelChemistry in ['B3LYP/cbsb7', 'B3LYP/6-311G(2d,d,p)', 'DFT_G03_b3lyp','B3LYP/6-311+G(3df,2p)']: + elif modelChemistry in ['B3LYP/cbsb7', 'B3LYP/6-311G(2d,d,p)', 'DFT_G03_b3lyp','B3LYP/6-311+G(3df,2p)','b3lyp/6-31G**']: bondEnergies = { 'C-H': 0.25, 'C-C': -1.89, 'C=C': -0.40, 'C#C': -1.50, 'O-H': -1.09, 'C-O': -1.18, 'C=O': -0.01, 'N-H': 1.36, 'C-N': -0.44, 'C#N': 0.22, 'C-S': -2.35, 'O=S': -5.19, 'S-H': -0.52, } @@ -727,6 +731,9 @@ def projectRotors(conformer, F, rotors, linear, TS): molecule `linear`, project out the nonvibrational modes from the force constant matrix and use this to determine the vibrational frequencies. The list of vibrational frequencies is returned in cm^-1. + + Refer to Gaussian whitepaper (http://gaussian.com/vib/) for procedure to calculate + harmonic oscillator vibrational frequencies using the force constant matrix. """ Nrotors = len(rotors) @@ -851,7 +858,9 @@ def projectRotors(conformer, F, rotors, linear, TS): logging.debug('Frequencies from internal Hessian') for i in range(3*Natoms-external): - logging.debug(numpy.sqrt(eig[i])/(2 * math.pi * constants.c * 100)) + with numpy.warnings.catch_warnings(): + numpy.warnings.filterwarnings('ignore', r'invalid value encountered in sqrt') + logging.debug(numpy.sqrt(eig[i])/(2 * math.pi * constants.c * 100)) # Now we can start thinking about projecting out the internal rotations Dint=numpy.zeros((3*Natoms,Nrotors), numpy.float64) @@ -930,7 +939,9 @@ def projectRotors(conformer, F, rotors, linear, TS): logging.debug('Frequencies from projected Hessian') for i in range(3*Natoms): - logging.debug(numpy.sqrt(eig[i])/(2 * math.pi * constants.c * 100)) + with numpy.warnings.catch_warnings(): + numpy.warnings.filterwarnings('ignore', r'invalid value encountered in sqrt') + logging.debug(numpy.sqrt(eig[i])/(2 * math.pi * constants.c * 100)) return numpy.sqrt(eig[-Nvib:]) / (2 * math.pi * constants.c * 100) diff --git a/rmgpy/cantherm/thermo.py b/rmgpy/cantherm/thermo.py index 1768322b10..1d8e00f3c0 100644 --- a/rmgpy/cantherm/thermo.py +++ b/rmgpy/cantherm/thermo.py @@ -164,7 +164,7 @@ def save(self, outputFile): f.write('{0}\n\n'.format(prettify(string))) f.close() - + # write chemkin file f = open(os.path.join(os.path.dirname(outputFile), 'chem.inp'), 'a') if isinstance(species, Species): if species.molecule and isinstance(species.molecule[0], Molecule): @@ -179,7 +179,14 @@ def save(self, outputFile): string = writeThermoEntry(species, elementCounts=elementCounts, verbose=False) f.write('{0}\n'.format(string)) f.close() - + + # write species dictionary + f = open(os.path.join(os.path.dirname(outputFile), 'species_dictionary.txt'), 'a') + if isinstance(species, Species): + if species.molecule and isinstance(species.molecule[0], Molecule): + f.write(species.molecule[0].toAdjacencyList(removeH=False,label=species.label)) + f.write('\n') + f.close() def plot(self, outputDirectory): """ diff --git a/scripts/sensitivity.py b/rmgpy/chemkin.pxd similarity index 60% rename from scripts/sensitivity.py rename to rmgpy/chemkin.pxd index 3d8421c8bb..f23407f709 100644 --- a/scripts/sensitivity.py +++ b/rmgpy/chemkin.pxd @@ -28,48 +28,7 @@ # # ############################################################################### -""" -This script runs stand-alone sensitivity analysis on an RMG job. -""" +from rmgpy.reaction cimport Reaction +from rmgpy.kinetics.model cimport KineticsModel -import os.path -import argparse - -from rmgpy.tools.sensitivity import runSensitivity - -################################################################################ - -def parse_arguments(): - - parser = argparse.ArgumentParser() - parser.add_argument('input', metavar='INPUT', type=str, nargs=1, - help='RMG input file') - parser.add_argument('chemkin', metavar='CHEMKIN', type=str, nargs=1, - help='Chemkin file') - parser.add_argument('dictionary', metavar='DICTIONARY', type=str, nargs=1, - help='RMG dictionary file') - parser.add_argument('--no-dlim', dest='dlim', action='store_false', - help='Turn off diffusion-limited rates for LiquidReactor') - args = parser.parse_args() - - inputFile = os.path.abspath(args.input[0]) - chemkinFile = os.path.abspath(args.chemkin[0]) - dictFile = os.path.abspath(args.dictionary[0]) - dflag = args.dlim - - return inputFile, chemkinFile, dictFile, dflag - -def main(): - # This might not work anymore because functions were modified for use with webserver - - inputFile, chemkinFile, dictFile, dflag = parse_arguments() - - runSensitivity(inputFile, chemkinFile, dictFile, dflag) - -################################################################################ - -if __name__ == '__main__': - main() - - - \ No newline at end of file +cpdef _process_duplicate_reactions(list reactionList) diff --git a/rmgpy/chemkin.py b/rmgpy/chemkin.pyx similarity index 97% rename from rmgpy/chemkin.py rename to rmgpy/chemkin.pyx index 7870fd2f08..6d470d0afe 100644 --- a/rmgpy/chemkin.py +++ b/rmgpy/chemkin.pyx @@ -841,13 +841,13 @@ def loadTransportFile(path, speciesDict): comment = comment.strip(), ) -def loadChemkinFile(path, dictionaryPath=None, transportPath=None, readComments = True, thermoPath = None, useChemkinNames=False): +def loadChemkinFile(path, dictionaryPath=None, transportPath=None, readComments=True, thermoPath=None, + useChemkinNames=False, checkDuplicates=True): """ Load a Chemkin input file located at `path` on disk to `path`, returning lists of the species and reactions in the Chemkin file. The 'thermoPath' point to a separate thermo file, or, if 'None' is specified, the function will look for the thermo database within the chemkin mechanism file """ - speciesList = []; speciesDict = {}; speciesAliases = {} reactionList = [] @@ -857,15 +857,13 @@ def loadChemkinFile(path, dictionaryPath=None, transportPath=None, readComments # HTML output. if dictionaryPath: speciesDict = loadSpeciesDictionary(dictionaryPath) - + with open(path, 'r+b') as f: line0 = f.readline() while line0 != '': line = removeCommentFromLine(line0)[0] line = line.strip() - tokens = line.split() - tokens_upper = line.upper().split() if 'SPECIES' in line.upper(): # Unread the line (we'll re-read it in readReactionBlock()) @@ -898,29 +896,69 @@ def loadChemkinFile(path, dictionaryPath=None, transportPath=None, readComments readThermoBlock(f, speciesDict) break line0 = f.readline() - # Index the reactions now to have identical numbering as in Chemkin + # Index the reactions now to have identical numbering as in Chemkin index = 0 for reaction in reactionList: index += 1 reaction.index = index - # Check for marked (and unmarked!) duplicate reactions - # Combine marked duplicate reactions into a single reaction using MultiKinetics - # Raise exception for unmarked duplicate reactions - duplicateReactionsToRemove = [] - duplicateReactionsToAdd = [] - for index1 in range(len(reactionList)): + # Process duplicate reactions + if checkDuplicates: + _process_duplicate_reactions(reactionList) + + # If the transport path is given, then read it to obtain the transport + # properties + if transportPath: + loadTransportFile(transportPath, speciesDict) + + if not useChemkinNames: + # Apply species aliases if known + for spec in speciesList: + try: + spec.label = speciesAliases[spec.label] + except KeyError: + pass + + # Attempt to extract index from species label + indexPattern = re.compile(r'\(\d+\)$') + for spec in speciesList: + if indexPattern.search(spec.label): + label, sep, index = spec.label[:-1].rpartition('(') + spec.label = label + spec.index = int(index) + + reactionList.sort(key=lambda reaction: reaction.index) + return speciesList, reactionList + + +cpdef _process_duplicate_reactions(list reactionList): + """ + Check for marked (and unmarked!) duplicate reactions + Combine marked duplicate reactions into a single reaction using MultiKinetics + Raise exception for unmarked duplicate reactions + """ + cdef list duplicateReactionsToRemove = [] + cdef list duplicateReactionsToAdd = [] + cdef int index1, index2 + cdef Reaction reaction, reaction1, reaction2 + cdef KineticsModel kinetics + + for index1 in xrange(len(reactionList)): reaction1 = reactionList[index1] if reaction1 in duplicateReactionsToRemove: continue - for index2 in range(index1+1, len(reactionList)): + for index2 in xrange(index1 + 1, len(reactionList)): reaction2 = reactionList[index2] - if reaction1.reactants == reaction2.reactants and reaction1.products == reaction2.products and reaction1.specificCollider == reaction2.specificCollider: + if (reaction1.reactants == reaction2.reactants + and reaction1.products == reaction2.products + and reaction1.specificCollider == reaction2.specificCollider): if reaction1.duplicate and reaction2.duplicate: - + if isinstance(reaction1, LibraryReaction) and isinstance(reaction2, LibraryReaction): - assert reaction1.library == reaction2.library + if reaction1.library != reaction2.library: + raise ChemkinError("Identical reactions {0} and {1} taken from different libraries: {2}, " + "{3}".format(reaction1, reaction2, reaction1.library, reaction2.library)) if reaction1 not in duplicateReactionsToRemove: # already created duplicate reaction, move on to appending any additional duplicate kinetics if isinstance(reaction1.kinetics, @@ -930,16 +968,19 @@ def loadChemkinFile(path, dictionaryPath=None, transportPath=None, readComments _kinetics.Arrhenius): kinetics = _kinetics.MultiArrhenius() else: - logging.warning('Unexpected kinetics type {0} for duplicate reaction {1}. Not combining reactions.'.format(reaction1.kinetics.__class__, reaction1)) + logging.warning( + 'Unexpected kinetics type {0} for duplicate reaction {1}. ' + 'Not combining reactions.'.format(reaction1.kinetics.__class__, reaction1) + ) continue reaction = LibraryReaction( - index = reaction1.index, - reactants = reaction1.reactants, - products = reaction1.products, - specificCollider = reaction1.specificCollider, - kinetics = kinetics, - library = reaction1.library, - duplicate = False, + index=reaction1.index, + reactants=reaction1.reactants, + products=reaction1.products, + specificCollider=reaction1.specificCollider, + kinetics=kinetics, + library=reaction1.library, + duplicate=False, ) duplicateReactionsToAdd.append(reaction) kinetics.arrhenius = [reaction1.kinetics] @@ -949,7 +990,7 @@ def loadChemkinFile(path, dictionaryPath=None, transportPath=None, readComments # Do not use as duplicate reactions if it's not a library reaction # Template reactions should be kept separate continue - + if (isinstance(reaction.kinetics, _kinetics.MultiPDepArrhenius) and isinstance(reaction2.kinetics, @@ -957,46 +998,22 @@ def loadChemkinFile(path, dictionaryPath=None, transportPath=None, readComments reaction.kinetics.arrhenius.append(reaction2.kinetics) elif (isinstance(reaction.kinetics, _kinetics.MultiArrhenius) and - isinstance(reaction2.kinetics, - _kinetics.Arrhenius)): + isinstance(reaction2.kinetics, + _kinetics.Arrhenius)): reaction.kinetics.arrhenius.append(reaction2.kinetics) else: raise ChemkinError('Mixed kinetics for duplicate reaction {0}.'.format(reaction)) - + duplicateReactionsToRemove.append(reaction2) elif reaction1.kinetics.isPressureDependent() == reaction2.kinetics.isPressureDependent(): - # If both reactions are pressure-independent or both are pressure-dependent, then they need duplicate tags - # Chemkin treates pdep and non-pdep reactions as different, so those are okay + # If both reactions are pressure-independent or both are pressure-dependent, then they need + # duplicate tags. Chemkin treates pdep and non-pdep reactions as different, so those are okay raise ChemkinError('Encountered unmarked duplicate reaction {0}.'.format(reaction1)) - + for reaction in duplicateReactionsToRemove: reactionList.remove(reaction) reactionList.extend(duplicateReactionsToAdd) - # If the transport path is given, then read it to obtain the transport - # properties - if transportPath: - loadTransportFile(transportPath, speciesDict) - - if not useChemkinNames: - # Apply species aliases if known - for spec in speciesList: - try: - spec.label = speciesAliases[spec.label] - except KeyError: - pass - - # Attempt to extract index from species label - indexPattern = re.compile(r'\(\d+\)$') - for spec in speciesList: - if indexPattern.search(spec.label): - label, sep, index = spec.label[:-1].rpartition('(') - spec.label = label - spec.index = int(index) - - reactionList.sort(key=lambda reaction: reaction.index) - return speciesList, reactionList - def readSpeciesBlock(f, speciesDict, speciesAliases, speciesList): """ diff --git a/rmgpy/chemkinTest.py b/rmgpy/chemkinTest.py index 790fc89215..b34b3a217b 100644 --- a/rmgpy/chemkinTest.py +++ b/rmgpy/chemkinTest.py @@ -32,11 +32,13 @@ import mock import os from chemkin import * -from chemkin import _removeLineBreaks +from chemkin import _removeLineBreaks, _process_duplicate_reactions import rmgpy from rmgpy.species import Species from rmgpy.reaction import Reaction -from rmgpy.kinetics.arrhenius import Arrhenius +from rmgpy.data.kinetics import LibraryReaction +from rmgpy.kinetics.arrhenius import Arrhenius, MultiArrhenius +from rmgpy.kinetics.chebyshev import Chebyshev ################################################### @@ -310,6 +312,76 @@ def testReadSpecificCollider(self): self.assertEqual(reaction.specificCollider.label, 'N2(5)') + def test_process_duplicate_reactions(self): + """ + Test that duplicate reactions are handled correctly when + loading a Chemkin file. + """ + s1 = Species().fromSMILES('CC') + s2 = Species().fromSMILES('[CH3]') + s3 = Species().fromSMILES('[OH]') + s4 = Species().fromSMILES('C[CH2]') + s5 = Species().fromSMILES('O') + r1 = Reaction(reactants=[s1], products=[s2, s2], duplicate=False, kinetics=Arrhenius()) + r2 = Reaction(reactants=[s1, s3], products=[s4, s5], duplicate=True, kinetics=Arrhenius()) + r3 = Reaction(reactants=[s1, s3], products=[s4, s5], duplicate=True, kinetics=Arrhenius()) + r4 = Reaction(reactants=[s1, s3], products=[s4, s5], duplicate=False, kinetics=Arrhenius()) + r5 = LibraryReaction(reactants=[s1, s3], products=[s4, s5], duplicate=True, + kinetics=Arrhenius(), library='lib1') + r6 = LibraryReaction(reactants=[s1, s3], products=[s4, s5], duplicate=True, + kinetics=Arrhenius(), library='lib2') + r7 = LibraryReaction(reactants=[s1, s3], products=[s4, s5], duplicate=True, + kinetics=Chebyshev(), library='lib1') + r8 = LibraryReaction(reactants=[s1, s3], products=[s4, s5], duplicate=True, + kinetics=Arrhenius(), library='lib1') + r9 = LibraryReaction(reactants=[s1, s3], products=[s4, s5], duplicate=False, + kinetics=MultiArrhenius(arrhenius=[Arrhenius(), Arrhenius()]), library='lib1') + reaction_list_with_duplicate = [r1, r2, r3] + reaction_list_with_duplicate2 = [r1, r2, r3] + reaction_list_unmarked_duplicate = [r1, r2, r4] + reaction_list_unequal_libraries = [r1, r5, r6] + reaction_list_mixed_kinetics = [r1, r5, r7] + reaction_list_mergeable = [r1, r5, r8] + reaction_list_merged = [r1, r9] + + # Test that duplicates are not removed for non-library reactions + _process_duplicate_reactions(reaction_list_with_duplicate) + self.assertEqual(reaction_list_with_duplicate, reaction_list_with_duplicate2) + + # Test that unmarked duplicate reactions are detected if both + # reactions are p-dep or p-indep + self.assertRaisesRegexp(ChemkinError, + 'Encountered unmarked duplicate reaction', + _process_duplicate_reactions, + reaction_list_unmarked_duplicate) + + # Test that unequal libraries are recognized + self.assertRaisesRegexp(ChemkinError, + 'from different libraries', + _process_duplicate_reactions, + reaction_list_unequal_libraries) + + # Test that an error is raised for reactions with kinetics + # that cannot be merged + self.assertRaisesRegexp(ChemkinError, + 'Mixed kinetics for duplicate reaction', + _process_duplicate_reactions, + reaction_list_mixed_kinetics) + + # Test that duplicate library reactions are merged successfully + _process_duplicate_reactions(reaction_list_mergeable) + self.assertEqual(len(reaction_list_mergeable), len(reaction_list_merged)) + self.assertEqual(reaction_list_mergeable[0], reaction_list_merged[0]) + rtest = reaction_list_mergeable[1] + rtrue = reaction_list_merged[1] + self.assertEqual(rtest.reactants, rtrue.reactants) + self.assertEqual(rtest.products, rtrue.products) + self.assertEqual(rtest.duplicate, rtrue.duplicate) + self.assertEqual(rtest.library, rtrue.library) + self.assertTrue(isinstance(rtest.kinetics, MultiArrhenius)) + self.assertTrue(all(isinstance(k, Arrhenius) for k in rtest.kinetics.arrhenius)) + + class TestReadReactionComments(unittest.TestCase): @classmethod def setUpClass(self): @@ -393,7 +465,7 @@ def setUpClass(self): Reaction index: Chemkin #2; RMG #4 Template reaction: R_Recombination Flux pairs: [CH3], CC; [CH3], CC; -From training reaction 21 for rate rule [C_rad/H2/Cs;C_methyl] +From training reaction 21 used for C_rad/H2/Cs;C_methyl Exact match found for rate rule [C_rad/H2/Cs;C_methyl] Euclidian distance = 0 """] diff --git a/rmgpy/data/base.py b/rmgpy/data/base.py index ff8fb96767..6b956873b7 100644 --- a/rmgpy/data/base.py +++ b/rmgpy/data/base.py @@ -1267,23 +1267,30 @@ def isMoleculeForbidden(self, molecule): contains forbidden functionality, or ``False`` if not. Labeled atoms on the forbidden structures and the molecule are honored. """ + from rmgpy.species import Species + for entry in self.entries.values(): - entryLabeledAtoms = entry.item.getLabeledAtoms() - moleculeLabeledAtoms = molecule.getLabeledAtoms() - initialMap = {} - for label in entryLabeledAtoms: - # all group labels must be present in the molecule - if label not in moleculeLabeledAtoms: break - initialMap[moleculeLabeledAtoms[label]] = entryLabeledAtoms[label] - else: - if molecule.isMappingValid(entry.item, initialMap) and molecule.isSubgraphIsomorphic(entry.item, initialMap): + if isinstance(entry.item, Molecule) or isinstance(entry.item, Species): + # Perform an isomorphism check + if entry.item.isIsomorphic(molecule): return True + elif isinstance(entry.item, Group): + # We need to do subgraph isomorphism + entryLabeledAtoms = entry.item.getLabeledAtoms() + moleculeLabeledAtoms = molecule.getLabeledAtoms() + initialMap = {} + for label in entryLabeledAtoms: + # all group labels must be present in the molecule + if label not in moleculeLabeledAtoms: break + initialMap[moleculeLabeledAtoms[label]] = entryLabeledAtoms[label] + else: + if molecule.isMappingValid(entry.item, initialMap) and molecule.isSubgraphIsomorphic(entry.item, initialMap): + return True + else: + raise NotImplementedError('Checking is only implemented for forbidden Groups, Molecule, and Species.') # Until we have more thermodynamic data of molecular ions we will forbid them - molecule_charge = 0 - for atom in molecule.atoms: - molecule_charge += atom.charge - if molecule_charge != 0: + if molecule.getNetCharge() != 0: return True return False @@ -1301,22 +1308,27 @@ def saveOld(self, path): """ self.saveOldDictionary(path) - def loadEntry(self, label, molecule=None, group=None, shortDesc='', longDesc=''): + def loadEntry(self, label, group=None, molecule=None, species=None, shortDesc='', longDesc=''): """ Load an entry from the forbidden structures database. This method is automatically called during loading of the forbidden structures database. """ - assert molecule is not None or group is not None - assert not (molecule is not None and group is not None) + from rmgpy.species import Species + + if sum([bool(molecule), bool(group), bool(species)]) != 1: + raise DatabaseError('A forbidden group should be defined with exactly one item from ' + 'the following options: group, molecule, or species.') if molecule is not None: - item = Molecule.fromAdjacencyList(molecule) + item = Molecule().fromAdjacencyList(molecule) + elif species is not None: + item = Species().fromAdjacencyList(species) + item.generate_resonance_structures() elif group is not None: - if ( group[0:3].upper() == 'OR{' or - group[0:4].upper() == 'AND{' or - group[0:7].upper() == 'NOT OR{' or - group[0:8].upper() == 'NOT AND{' - ): + if (group[0:3].upper() == 'OR{' or + group[0:4].upper() == 'AND{' or + group[0:7].upper() == 'NOT OR{' or + group[0:8].upper() == 'NOT AND{'): item = makeLogicNode(group) else: item = Group().fromAdjacencyList(group) diff --git a/rmgpy/data/baseTest.py b/rmgpy/data/baseTest.py index c1f4528c6e..9fab6afa90 100644 --- a/rmgpy/data/baseTest.py +++ b/rmgpy/data/baseTest.py @@ -32,8 +32,8 @@ import unittest from external.wip import work_in_progress -from rmgpy.data.base import Entry, Database -from rmgpy.molecule import Group +from rmgpy.data.base import Entry, Database, ForbiddenStructures +from rmgpy.molecule import Group, Molecule ################################################################################ @@ -117,6 +117,103 @@ def testMatchNodeToNode(self): self.assertFalse(self.database.matchNodeToNode(entry1,entry2)) +class TestForbiddenStructures(unittest.TestCase): + + def setUp(self): + self.database = ForbiddenStructures() + + def test_forbidden_group(self): + """Test that we can load and check a forbidden group.""" + test = """ +1 C u2 p0 {2,D} +2 C u0 {1,D} +""" + self.database.loadEntry( + label='test', + group=test, + ) + + molecule = Molecule().fromAdjacencyList(""" +multiplicity 3 +1 C u2 p0 c0 {2,D} +2 C u0 p0 c0 {1,D} {3,S} {4,S} +3 H u0 p0 c0 {2,S} +4 H u0 p0 c0 {2,S} +""") + + self.assertTrue(self.database.isMoleculeForbidden(molecule)) + + def test_forbidden_molecule(self): + """Test that we can load and check a forbidden molecule.""" + test = """ +1 C u4 p0 c0 +""" + self.database.loadEntry( + label='test', + molecule=test, + ) + + molecule = Molecule().fromAdjacencyList(""" +1 C u4 p0 c0 +""") + + self.assertTrue(self.database.isMoleculeForbidden(molecule)) + + def test_forbidden_species(self): + """Test that we can load and check a forbidden species. + + This is a hypothetical test, we don't actually forbid benzene.""" + test = """ +1 C u0 p0 c0 {2,D} {6,S} {7,S} +2 C u0 p0 c0 {1,D} {3,S} {8,S} +3 C u0 p0 c0 {2,S} {4,D} {9,S} +4 C u0 p0 c0 {3,D} {5,S} {10,S} +5 C u0 p0 c0 {4,S} {6,D} {11,S} +6 C u0 p0 c0 {1,S} {5,D} {12,S} +7 H u0 p0 c0 {1,S} +8 H u0 p0 c0 {2,S} +9 H u0 p0 c0 {3,S} +10 H u0 p0 c0 {4,S} +11 H u0 p0 c0 {5,S} +12 H u0 p0 c0 {6,S} +""" + self.database.loadEntry( + label='test', + species=test, + ) + + molecule1 = Molecule().fromAdjacencyList(""" +1 C u0 p0 c0 {2,D} {6,S} {7,S} +2 C u0 p0 c0 {1,D} {3,S} {8,S} +3 C u0 p0 c0 {2,S} {4,D} {9,S} +4 C u0 p0 c0 {3,D} {5,S} {10,S} +5 C u0 p0 c0 {4,S} {6,D} {11,S} +6 C u0 p0 c0 {1,S} {5,D} {12,S} +7 H u0 p0 c0 {1,S} +8 H u0 p0 c0 {2,S} +9 H u0 p0 c0 {3,S} +10 H u0 p0 c0 {4,S} +11 H u0 p0 c0 {5,S} +12 H u0 p0 c0 {6,S} +""") + molecule2 = Molecule().fromAdjacencyList(""" +1 C u0 p0 c0 {2,B} {6,B} {7,S} +2 C u0 p0 c0 {1,B} {3,B} {8,S} +3 C u0 p0 c0 {2,B} {4,B} {9,S} +4 C u0 p0 c0 {3,B} {5,B} {10,S} +5 C u0 p0 c0 {4,B} {6,B} {11,S} +6 C u0 p0 c0 {1,B} {5,B} {12,S} +7 H u0 p0 c0 {1,S} +8 H u0 p0 c0 {2,S} +9 H u0 p0 c0 {3,S} +10 H u0 p0 c0 {4,S} +11 H u0 p0 c0 {5,S} +12 H u0 p0 c0 {6,S} +""") + + self.assertTrue(self.database.isMoleculeForbidden(molecule1)) + self.assertTrue(self.database.isMoleculeForbidden(molecule2)) + ################################################################################ diff --git a/rmgpy/data/kinetics/common.py b/rmgpy/data/kinetics/common.py index a6a320acce..314e8d6adc 100644 --- a/rmgpy/data/kinetics/common.py +++ b/rmgpy/data/kinetics/common.py @@ -156,8 +156,8 @@ def filter_reactions(reactants, products, reactionList): warnings.warn("The filter_reactions method is no longer used and may be removed in a future version.", DeprecationWarning) # Convert from molecules to species and generate resonance isomers. - reactants = ensure_species(reactants, resonance=True) - products = ensure_species(products, resonance=True) + ensure_species(reactants, resonance=True) + ensure_species(products, resonance=True) reactions = reactionList[:] @@ -197,11 +197,10 @@ def filter_reactions(reactants, products, reactionList): def ensure_species(input_list, resonance=False, keepIsomorphic=False): """ - Given an input list of molecules or species, return a list with only - species objects. + The input list of :class:`Species` or :class:`Molecule` objects is modified + in place to only have :class:`Species` objects. Returns None. """ - output_list = [] - for item in input_list: + for index, item in enumerate(input_list): if isinstance(item, Molecule): new_item = Species(molecule=[item]) elif isinstance(item, Species): @@ -210,9 +209,7 @@ def ensure_species(input_list, resonance=False, keepIsomorphic=False): raise TypeError('Only Molecule or Species objects can be handled.') if resonance: new_item.generate_resonance_structures(keepIsomorphic=keepIsomorphic) - output_list.append(new_item) - - return output_list + input_list[index] = new_item def generate_molecule_combos(input_species): @@ -231,13 +228,15 @@ def generate_molecule_combos(input_species): def ensure_independent_atom_ids(input_species, resonance=True): """ - Given a list or tuple of :class:`Species` objects, ensure that atom ids are - independent across all of the species. Optionally, the `resonance` argument - can be set to False to not generate resonance structures. + Given a list or tuple of :class:`Species` or :class:`Molecule` objects, + ensure that atom ids are independent. + The `resonance` argument can be set to False to not generate + resonance structures. - Modifies the input species in place, nothing is returned. + Modifies the list in place (replacing :class:`Molecule` with :class:`Species`). + Returns None. """ - + ensure_species(input_species, resonance=resonance) # Method to check that all species' atom ids are different def independent_ids(): num_atoms = 0 diff --git a/rmgpy/data/kinetics/database.py b/rmgpy/data/kinetics/database.py index 77f8778116..c0303d2d09 100644 --- a/rmgpy/data/kinetics/database.py +++ b/rmgpy/data/kinetics/database.py @@ -364,7 +364,7 @@ def generate_reactions(self, reactants, products=None, only_families=None, reson reactionList = [] if only_families is None: reactionList.extend(self.generate_reactions_from_libraries(reactants, products)) - reactionList.extend(self.generate_reactions_from_families(reactants, products, only_families=None, resonance=True)) + reactionList.extend(self.generate_reactions_from_families(reactants, products, only_families=None, resonance=resonance)) return reactionList def generate_reactions_from_libraries(self, reactants, products=None): @@ -386,7 +386,7 @@ def generate_reactions_from_library(self, library, reactants, products=None): provided `reactants`, which can be either :class:`Molecule` objects or :class:`Species` objects. """ - reactants = ensure_species(reactants) + ensure_species(reactants) reaction_list = [] for entry in library.entries.values(): @@ -434,10 +434,9 @@ def generate_reactions_from_families(self, reactants, products=None, only_famili elif reactants[0].isIsomorphic(reactants[1]): same_reactants = True - # Convert to Species objects if necessary - reactants = ensure_species(reactants) - - # Label reactant atoms for proper degeneracy calculation + # Label reactant atoms for proper degeneracy calculation (cannot be in tuple) + if isinstance(reactants, tuple): + reactants = list(reactants) ensure_independent_atom_ids(reactants, resonance=resonance) combos = generate_molecule_combos(reactants) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 09fbc1d173..d839063abd 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -46,7 +46,8 @@ from rmgpy.molecule.resonance import generate_aromatic_resonance_structures from rmgpy.species import Species -from .common import saveEntry, ensure_species, find_degenerate_reactions, generate_molecule_combos +from .common import saveEntry, ensure_species, find_degenerate_reactions, generate_molecule_combos,\ + ensure_independent_atom_ids from .depository import KineticsDepository from .groups import KineticsGroups from .rules import KineticsRules @@ -1067,7 +1068,7 @@ def addKineticsRulesFromTrainingSet(self, thermoDatabase=None): shortDesc="Rate rule generated from training reaction {0}. ".format(entry.index) + entry.shortDesc, longDesc="Rate rule generated from training reaction {0}. ".format(entry.index) + entry.longDesc, ) - new_entry.data.comment = "From training reaction {1} for rate rule {0}".format(';'.join([g.label for g in template]), entry.index) + new_entry.data.comment = "From training reaction {1} used for {0}".format(';'.join([g.label for g in template]), entry.index) new_entry.data.A.value_si /= entry.item.degeneracy try: @@ -1115,7 +1116,7 @@ def addKineticsRulesFromTrainingSet(self, thermoDatabase=None): shortDesc="Rate rule generated from training reaction {0}. ".format(entry.index) + entry.shortDesc, longDesc="Rate rule generated from training reaction {0}. ".format(entry.index) + entry.longDesc, ) - new_entry.data.comment = "From training reaction {1} for rate rule {0}".format(';'.join([g.label for g in template]), entry.index) + new_entry.data.comment = "From training reaction {1} used for {0}".format(';'.join([g.label for g in template]), entry.index) new_entry.data.A.value_si /= new_degeneracy try: @@ -1451,14 +1452,14 @@ def generateReactions(self, reactants, products=None, prod_resonance=True): consistent with the template of this reaction family. Args: - reactants: List of Molecules to react - products: List of Molecules or Species of desired product structures (optional) - prod_resonance: Flag to generate resonance structures for product checking (optional) - Defaults to True, resonance structures are compared + reactants (list): List of Molecules to react. + products (list, optional): List of Molecules or Species of desired product structures. + prod_resonance (bool, optional): Flag to generate resonance structures for product checking. + Defaults to True, resonance structures are compared. Returns: List of all reactions containing Molecule objects with the - specified reactants and products within this family. + specified reactants and products within this family. Degenerate reactions are returned as separate reactions. """ reactionList = [] @@ -1545,22 +1546,28 @@ def calculateDegeneracy(self, reaction): `ignoreSameReactants= True` to this method. """ reaction.degeneracy = 1 - - # find combinations of resonance isomers - specReactants = ensure_species(reaction.reactants, resonance=True, keepIsomorphic=True) - molecule_combos = generate_molecule_combos(specReactants) + # Check if the reactants are the same + # If they refer to the same memory address, then make a deep copy so + # they can be manipulated independently + reactants = reaction.reactants + same_reactants = False + if len(reactants) == 2: + if reactants[0] is reactants[1]: + reactants[1] = reactants[1].copy(deep=True) + same_reactants = True + elif reactants[0].isIsomorphic(reactants[1]): + same_reactants = True + + # Label reactant atoms for proper degeneracy calculation + ensure_independent_atom_ids(reactants, resonance=True) + molecule_combos = generate_molecule_combos(reactants) reactions = [] for combo in molecule_combos: reactions.extend(self.__generateReactions(combo, products=reaction.products, forward=True)) - # Check if the reactants are the same - sameReactants = False - if len(specReactants) == 2 and specReactants[0].isIsomorphic(specReactants[1]): - sameReactants = True - # remove degenerate reactions - reactions = find_degenerate_reactions(reactions, sameReactants, kinetics_family=self) + reactions = find_degenerate_reactions(reactions, same_reactants, kinetics_family=self) # remove reactions with different templates (only for TemplateReaction) if isinstance(reaction, TemplateReaction): @@ -1692,7 +1699,7 @@ def __generateReactions(self, reactants, products=None, forward=True, prod_reson # If products is given, remove reactions from the reaction list that # don't generate the given products if products is not None: - products = ensure_species(products, resonance=prod_resonance) + ensure_species(products, resonance=prod_resonance) rxnList0 = rxnList[:] rxnList = [] @@ -1836,7 +1843,8 @@ def getReactionPairs(self, reaction): error = True if error: - raise ReactionPairsError('Unable to determine reaction pairs for {0!s} reaction {1!s}.'.format(self.label, reaction)) + logging.debug('Preset mapping missing for determining reaction pairs for family {0!s}, falling back to Reaction.generatePairs'.format(self.label)) + return [] else: return pairs @@ -1913,9 +1921,9 @@ def getKinetics(self, reaction, templateLabels, degeneracy=1, estimator='', retu 2. the source - this will be `None` if from a template estimate 3. the entry - this will be `None` if from a template estimate 4. isForward a boolean denoting whether the matched entry is in the same - direction as the inputted reaction. This will always be True if using - rates rules or group additivity. This can be `True` or `False` if using - a depository + direction as the inputted reaction. This will always be True if using + rates rules or group additivity. This can be `True` or `False` if using + a depository If returnAllKinetics==False, only the first (best?) matching kinetics is returned. """ @@ -2124,6 +2132,15 @@ def addAtomLabelsForReaction(self, reaction, output_with_resonance = True): # make sure we start with reaction with species objects reaction.ensure_species(reactant_resonance=False, product_resonance=False) + reactants = reaction.reactants + products = reaction.products + # ensure all species are independent references + if len(reactants + products) > len(set([id(s) for s in reactants + products])): + logging.debug('Copying reactants and products for reaction {} since they have identical species references'.format(reaction)) + # not all species are independent + reactants = [s.copy(deep=True) for s in reactants] + products = [s.copy(deep=True) for s in products] + # get all possible pairs of resonance structures reactant_pairs = list(itertools.product(*[s.molecule for s in reaction.reactants])) product_pairs = list(itertools.product(*[s.molecule for s in reaction.products])) @@ -2144,20 +2161,22 @@ def addAtomLabelsForReaction(self, reaction, output_with_resonance = True): # place the molecules in reaction's species object # this prevents overwriting of attributes of species objects by this method - for species in reaction.products: + for index, species in enumerate(products): for labeled_molecule in labeled_products: if species.isIsomorphic(labeled_molecule): species.molecule = [labeled_molecule] + reaction.products[index] = species break else: - raise ActionError('Could not find isomorphic molecule to fit the original product {}'.format(species)) - for species in reaction.reactants: + raise ActionError('Could not find isomorphic molecule to fit the original product {} from reaction {}'.format(species, reaction)) + for index, species in enumerate(reactants): for labeled_molecule in labeled_reactants: if species.isIsomorphic(labeled_molecule): species.molecule = [labeled_molecule] + reaction.reactants[index] = species break else: - raise ActionError('Could not find isomorphic molecule to fit the original reactant {}'.format(species)) + raise ActionError('Could not find isomorphic molecule to fit the original reactant {} from reaction {}'.format(species, reaction)) if output_with_resonance: # convert the molecules to species objects with resonance structures @@ -2299,21 +2318,22 @@ def extractSourceFromComments(self, reaction): Will return the template associated with the matched rate rule. Returns a tuple containing (Boolean_Is_Kinetics_From_Training_reaction, Source_Data) - For a training reaction, the Source_Data returns - [Family_Label, Training_Reaction_Entry, Kinetics_In_Reverse?] - - For a reaction from rate rules, the Source_Data is a tuple containing - [Family_Label, {'template': originalTemplate, - 'degeneracy': degeneracy, - 'exact': boolean_exact?, - 'rules': a list of (original rate rule entry, weight in average) - 'training': a list of (original rate rule entry associated with training entry, original training entry, weight in average) - }] - where TrainingReactions are ones that have created rules used in the estimate. + For a training reaction, the Source_Data returns:: + + [Family_Label, Training_Reaction_Entry, Kinetics_In_Reverse?] - where Exact is a boolean of whether the rate is an exact match, - Template is the reaction template used, - and RateRules is a list of the rate rule entries containing the kinetics used + For a reaction from rate rules, the Source_Data is a tuple containing:: + + [Family_Label, {'template': originalTemplate, + 'degeneracy': degeneracy, + 'exact': boolean_exact?, + 'rules': a list of (original rate rule entry, weight in average) + 'training': a list of (original rate rule entry associated with training entry, original training entry, weight in average)}] + + + where Exact is a boolean of whether the rate is an exact match, Template is + the reaction template used, RateRules is a list of the rate rule entries containing + the kinetics used, and TrainingReactions are ones that have created rules used in the estimate. """ import re lines = reaction.kinetics.comment.split('\n') @@ -2354,7 +2374,12 @@ def extractSourceFromComments(self, reaction): # The rate rule string is right after the phrase 'for rate rule' rateRuleString = fullCommentString.split("for rate rule",1)[1].split()[0] - templateLabel = re.split(regex, rateRuleString)[1] + + if rateRuleString[0] == '[': + templateLabel = re.split(regex, rateRuleString)[1] + else: + templateLabel = rateRuleString #if has the line 'From training reaction # for rate rule node1;node2' + template = self.retrieveTemplate(templateLabel.split(';')) rules, trainingEntries = self.getSourcesForTemplate(template) diff --git a/rmgpy/data/kinetics/groups.py b/rmgpy/data/kinetics/groups.py index 8c051e7e77..a42a3658ea 100644 --- a/rmgpy/data/kinetics/groups.py +++ b/rmgpy/data/kinetics/groups.py @@ -154,16 +154,24 @@ def getReactionTemplate(self, reaction): # template is a list of the actual matched nodes # forwardTemplate is a list of the top level nodes that should be matched if len(template) != len(forwardTemplate): -# print 'len(template):', len(template) -# print 'len(forwardTemplate):', len(forwardTemplate) - msg = 'Unable to find matching template for reaction {0} in reaction family {1}.'.format(str(reaction), str(self)) + msg = 'Unable to find matching template for reaction {0} in reaction family {1}.'.format(str(reaction), str(self)) msg += 'Trying to match {0} but matched {1}'.format(str(forwardTemplate),str(template)) -# print 'reactants' -# for reactant in reaction.reactants: -# print reactant.toAdjacencyList() + '\n' -# print 'products' -# for product in reaction.products: -# print product.toAdjacencyList() + '\n' + logging.debug('len(template): {0}'.format(len(template))) + logging.debug('len(forwardTemplate): {0}'.format(len(forwardTemplate))) + logging.debug('reactants:') + for reactant in reaction.reactants: + if isinstance(reactant,Species): + for mol in reactant.molecule: + logging.debug(mol.toAdjacencyList()) + else: + logging.debug(reactant.toAdjacencyList()) + logging.debug('products:') + for product in reaction.products: + if isinstance(product,Species): + for mol in product.molecule: + logging.debug(mol.toAdjacencyList()) + else: + logging.debug(product.toAdjacencyList()) raise UndeterminableKineticsError(reaction, message=msg) return template diff --git a/rmgpy/data/kinetics/kineticsTest.py b/rmgpy/data/kinetics/kineticsTest.py index c7b9ebef73..5c09f4efbd 100644 --- a/rmgpy/data/kinetics/kineticsTest.py +++ b/rmgpy/data/kinetics/kineticsTest.py @@ -688,6 +688,21 @@ def test_ensure_independent_atom_ids(self): # checks second resonance structure id self.assertNotEqual(s2.molecule[1].atoms[0].id, -1) + def test_ensure_independent_atom_ids_no_resonance(self): + """ + Ensure ensure_independent_atom_ids does not generate resonance + """ + s1 = Species().fromSMILES('CCC') + s2 = Species().fromSMILES('C=C[CH]C') + self.assertEqual(s2.molecule[0].atoms[0].id, -1) + + ensure_independent_atom_ids([s1, s2],resonance=False) + # checks resonance structures + self.assertEqual(len(s2.molecule),1) + # checks that atom ids are changed + for atom in s2.molecule[0].atoms: + self.assertNotEqual(atom.id, -1) + def testSaveEntry(self): """ tests that save entry can run @@ -839,25 +854,171 @@ def test_generate_reactions_from_families_product_resonance(self): self.assertEqual(len(reaction_list), 1) self.assertEqual(reaction_list[0].degeneracy, 2) - reaction_list = self.database.kinetics.generate_reactions_from_families(reactants, products, only_families=['H_Abstraction'], resonance=False) + + def test_generate_reactions_from_families_product_resonance2(self): + """Test that we can specify the no product resonance structure when generating reactions""" + reactants = [ + Molecule().fromSMILES('CCC=C'), + Molecule().fromSMILES('[H]'), + ] + products = [ + Molecule().fromSMILES('CC=C[CH2]'), + Molecule().fromSMILES('[H][H]'), + ] + + reaction_list = self.database.kinetics.generate_reactions_from_families(reactants, products, only_families=['H_Abstraction'], resonance=False) self.assertEqual(len(reaction_list), 0) + self.assertTrue(isinstance(products[0],Species)) + self.assertEqual(len(products[0].molecule),1) + def test_generate_reactions_from_libraries(self): """Test that we can generate reactions from libraries""" reactants = [ Molecule().fromSMILES('CC=O'), Molecule().fromSMILES('[H]'), ] - products = [ - Molecule().fromSMILES('[CH2]C=O'), - Molecule().fromSMILES('[H][H]'), - ] reaction_list = self.database.kinetics.generate_reactions_from_libraries(reactants) self.assertEqual(len(reaction_list), 3) + def test_generate_reactions_from_libraries2(self): + """Test that we can generate reactions from libraries specifying products""" + reactants = [ + Molecule().fromSMILES('CC=O'), + Molecule().fromSMILES('[H]'), + ] + products = [ + Molecule().fromSMILES('[CH2]C=O'), + Molecule().fromSMILES('[H][H]'), + ] reaction_list_2 = self.database.kinetics.generate_reactions_from_libraries(reactants, products) self.assertEqual(len(reaction_list_2), 1) + + def test_add_atom_labels_for_reaction(self): + """Test that addAtomLabelsForReaction can identify reactions with resonance + The molecule [CH]=C=C has resonance in this reaction""" + from rmgpy.data.rmg import getDB + reactants = [ + Molecule().fromSMILES('C=C=C'), + Molecule().fromSMILES('[CH]=C=C'), + ] + products = [ + Molecule().fromSMILES('C#C[CH2]'), + Molecule().fromSMILES('C#CC'), + ] + reaction = TemplateReaction(reactants =reactants, + products = products, + family = 'H_Abstraction') + reaction.ensure_species(reactant_resonance=True, product_resonance=True) + family = getDB('kinetics').families['H_Abstraction'] + family.addAtomLabelsForReaction(reaction, output_with_resonance=False) + + # test that the reaction has labels + found_labels = [] + for species in reaction.reactants: + for atom in species.molecule[0].atoms: + if atom.label != '': + found_labels.append(atom.label) + self.assertEqual(len(found_labels), 3) + self.assertIn('*1',found_labels) + self.assertIn('*2',found_labels) + self.assertIn('*3',found_labels) + + # test for the products too + found_labels = [] + for species in reaction.products: + for atom in species.molecule[0].atoms: + if atom.label != '': + found_labels.append(atom.label) + self.assertEqual(len(found_labels), 3) + self.assertIn('*1',found_labels) + self.assertIn('*2',found_labels) + self.assertIn('*3',found_labels) + + def test_add_atom_labels_for_reaction_2(self): + """Test that addAtomLabelsForReaction can identify reactions with identical references + The molecule [CH]=C=C has resonance in this reaction""" + from rmgpy.data.rmg import getDB + s1 = Species().fromSMILES('C=C=C') + s2 = Species().fromSMILES('C=C=[CH]') + s3 = Species().fromSMILES('C#CC') + s2.generate_resonance_structures() + reactants = [s1,s2] + products = [s2,s3] + reaction = TemplateReaction(reactants =reactants, + products = products, + family = 'H_Abstraction') + family = getDB('kinetics').families['H_Abstraction'] + print reaction.reactants + print reaction.products + family.addAtomLabelsForReaction(reaction, output_with_resonance=False) + + # test that the reaction has labels + found_labels = [] + for species in reaction.reactants: + for atom in species.molecule[0].atoms: + if atom.label != '': + found_labels.append(atom.label) + self.assertEqual(len(found_labels), 3,'wrong number of labels found {0}'.format(found_labels)) + self.assertIn('*1',found_labels) + self.assertIn('*2',found_labels) + self.assertIn('*3',found_labels) + + # test for the products too + found_labels = [] + for species in reaction.products: + for atom in species.molecule[0].atoms: + if atom.label != '': + found_labels.append(atom.label) + self.assertEqual(len(found_labels), 3) + self.assertIn('*1',found_labels) + self.assertIn('*2',found_labels) + self.assertIn('*3',found_labels) + + def test_add_atom_labels_for_reaction_3(self): + """Test that addAtomLabelsForReaction can identify reactions with resonance and isotopes""" + from rmgpy.data.rmg import getDB + mr0 = Molecule().fromAdjacencyList('1 C u0 p0 c0 i13 {3,D} {4,S} {5,S}\n2 *1 C u0 p0 c0 {3,D} {6,S} {7,S}\n3 C u0 p0 c0 {1,D} {2,D}\n4 H u0 p0 c0 {1,S}\n5 H u0 p0 c0 {1,S}\n6 H u0 p0 c0 {2,S}\n7 *4 H u0 p0 c0 {2,S}\n') + mr1a = Molecule().fromAdjacencyList('multiplicity 2\n1 C u0 p0 c0 i13 {2,D} {4,S} {5,S}\n2 C u0 p0 c0 {1,D} {3,D}\n3 *1 C u1 p0 c0 {2,D} {6,S}\n4 H u0 p0 c0 {1,S}\n5 H u0 p0 c0 {1,S}\n6 H u0 p0 c0 {3,S}\n') + mr1b = Molecule().fromAdjacencyList('multiplicity 2\n1 C u1 p0 c0 i13 {2,S} {4,S} {5,S}\n2 C u0 p0 c0 {1,S} {3,T}\n3 *1 C u0 p0 c0 {2,T} {6,S}\n4 H u0 p0 c0 {1,S}\n5 H u0 p0 c0 {1,S}\n6 H u0 p0 c0 {3,S}\n') + mp1a = Molecule().fromAdjacencyList('multiplicity 2\n1 C u0 p0 c0 {2,D} {4,S} {5,S}\n2 C u0 p0 c0 {1,D} {3,D}\n3 *1 C u1 p0 c0 i13 {2,D} {6,S}\n4 H u0 p0 c0 {1,S}\n5 H u0 p0 c0 {1,S}\n6 H u0 p0 c0 {3,S}\n') + mp1b = Molecule().fromAdjacencyList('multiplicity 2\n1 C u1 p0 c0 {2,S} {4,S} {5,S}\n2 C u0 p0 c0 {1,S} {3,T}\n3 *1 C u0 p0 c0 i13 {2,T} {6,S}\n4 H u0 p0 c0 {1,S}\n5 H u0 p0 c0 {1,S}\n6 H u0 p0 c0 {3,S}\n') + s1 = Species(molecule = [mr0]) + s2 = Species(molecule = [mr1a,mr1b]) + s3 = Species(molecule = [mp1a,mp1b]) + reactants = [s1,s2] + products = [s1,s3] + reaction = TemplateReaction(reactants =reactants, + products = products, + family = 'H_Abstraction') + family = getDB('kinetics').families['H_Abstraction'] + print reaction.reactants + print reaction.products + family.addAtomLabelsForReaction(reaction, output_with_resonance=False) + + # test that the reaction has labels + found_labels = [] + for species in reaction.reactants: + for atom in species.molecule[0].atoms: + if atom.label != '': + found_labels.append(atom.label) + self.assertEqual(len(found_labels), 3,'wrong number of labels found {0}'.format(found_labels)) + self.assertIn('*1',found_labels) + self.assertIn('*2',found_labels) + self.assertIn('*3',found_labels) + + # test for the products too + found_labels = [] + for species in reaction.products: + for atom in species.molecule[0].atoms: + if atom.label != '': + found_labels.append(atom.label) + self.assertEqual(len(found_labels), 3) + self.assertIn('*1',found_labels) + self.assertIn('*2',found_labels) + self.assertIn('*3',found_labels) + diff --git a/rmgpy/data/kinetics/library.py b/rmgpy/data/kinetics/library.py index a71395ecf6..24b8ba59b8 100644 --- a/rmgpy/data/kinetics/library.py +++ b/rmgpy/data/kinetics/library.py @@ -152,10 +152,9 @@ def getLibraryReactions(self): c = entry._longDesc.split('\n') family_comments = [i for i in c if 'family: ' in i] familyname = family_comments[0].replace('family: ','') - logging.info(familyname) tstring = c[0] ind = tstring.find('rate rule') - tstring = tstring[ind+9:] + tstring = tstring[ind+10:] tstrings = tstring.split(';') tstrings[0] = tstrings[0][1:] tstrings[-1] = tstrings[-1][:-1] diff --git a/rmgpy/data/thermo.py b/rmgpy/data/thermo.py index 5516c37f35..9e718b21b8 100644 --- a/rmgpy/data/thermo.py +++ b/rmgpy/data/thermo.py @@ -108,6 +108,9 @@ def saveEntry(f, entry): f.write(' ],\n') if entry.data.Tmin is not None: f.write(' Tmin = {0!r},\n'.format(entry.data.Tmin)) if entry.data.Tmax is not None: f.write(' Tmax = {0!r},\n'.format(entry.data.Tmax)) + if entry.data.E0 is not None: f.write(' E0 = {0!r},\n'.format(entry.data.E0)) + if entry.data.Cp0 is not None: f.write(' Cp0 = {0!r},\n'.format(entry.data.Cp0)) + if entry.data.CpInf is not None: f.write(' CpInf = {0!r},\n'.format(entry.data.CpInf)) f.write(' ),\n') else: f.write(' thermo = {0!r},\n'.format(entry.data)) @@ -1292,10 +1295,15 @@ def getAllThermoData(self, species): thermoDataList.append(data) # Last entry is always the estimate from group additivity # Make it a tuple - data = (self.getThermoDataFromGroups(species), None, None) - # update group activity for symmetry - data[0].S298.value_si -= constants.R * math.log(species.getSymmetryNumber()) - thermoDataList.append(data) + try: + data = (self.getThermoDataFromGroups(species), None, None) + except DatabaseError: + # We don't have a GAV estimate, e.g. unsupported element + pass + else: + # update group activity for symmetry + data[0].S298.value_si -= constants.R * math.log(species.getSymmetryNumber()) + thermoDataList.append(data) # Return all of the resulting thermo parameters return thermoDataList diff --git a/rmgpy/data/thermoTest.py b/rmgpy/data/thermoTest.py index 2ac3802892..0f6aa4e812 100644 --- a/rmgpy/data/thermoTest.py +++ b/rmgpy/data/thermoTest.py @@ -314,6 +314,18 @@ def testThermoEstimationNotAffectDatabase(self): self.assertAlmostEqual(previous_enthalpy, latter_enthalpy, 2) + def test_getAllThermoData_fails_quietly(self): + """Test that getAllThermoData doesn't break when GAV fails.""" + spec = Species().fromSMILES('[Ne]') + + # Check that GAV fails + with self.assertRaises(DatabaseError): + self.database.getThermoDataFromGroups(spec) + + # Check that getAllThermoData doesn't break + thermo = self.database.getAllThermoData(spec) + self.assertEqual(len(thermo), 1) + class TestThermoAccuracy(unittest.TestCase): """ diff --git a/rmgpy/exceptions.py b/rmgpy/exceptions.py index ee93662c2e..27df3cd937 100644 --- a/rmgpy/exceptions.py +++ b/rmgpy/exceptions.py @@ -268,6 +268,12 @@ class VF2Error(Exception): """ pass +class CoreError(Exception): + """ + An exception raised if there is a problem within the model core + """ + pass + ################## move classes that extend off previous exceptions here class InvalidMicrocanonicalRateError(NetworkError): diff --git a/rmgpy/kinetics/arrhenius.pyx b/rmgpy/kinetics/arrhenius.pyx index 3bd7d811c6..1fd764af83 100644 --- a/rmgpy/kinetics/arrhenius.pyx +++ b/rmgpy/kinetics/arrhenius.pyx @@ -32,6 +32,7 @@ from libc.math cimport exp, log, sqrt, log10 cimport rmgpy.constants as constants import rmgpy.quantity as quantity +from rmgpy.exceptions import KineticsError ################################################################################ cdef class Arrhenius(KineticsModel): @@ -140,6 +141,10 @@ cdef class Arrhenius(KineticsModel): """ import numpy.linalg import scipy.stats + + assert len(Tlist) == len(klist), "length of temperatures and rates must be the same" + if len(Tlist) < 3+threeParams: + raise KineticsError('Not enough degrees of freedom to fit this Arrhenius expression') if threeParams: A = numpy.zeros((len(Tlist),3), numpy.float64) A[:,0] = numpy.ones_like(Tlist) diff --git a/rmgpy/kinetics/chebyshev.pyx b/rmgpy/kinetics/chebyshev.pyx index 1e45c93c33..7e598eb7f2 100644 --- a/rmgpy/kinetics/chebyshev.pyx +++ b/rmgpy/kinetics/chebyshev.pyx @@ -32,7 +32,8 @@ from libc.math cimport exp, log, sqrt, log10 cimport rmgpy.constants as constants import rmgpy.quantity as quantity - +import logging +from rmgpy.exceptions import KineticsError ################################################################################ cdef class Chebyshev(PDepKineticsModel): @@ -181,6 +182,18 @@ cdef class Chebyshev(PDepKineticsModel): cdef int t1, p1, t2, p2 cdef double T, P + if nT <= degreeT or nP <= degreeP: + raise KineticsError( + "The master equation data needs more temperature and pressure data "\ + "points than are placed into Chebyshev polynomial. Currently, the "\ + "data has {0} temperatures and the polynomial is set to have {1}. "\ + "The data has {2} pressures and the polynomial is set ot have {3}"\ + "".format(nT,degreeT,nP,degreeP)) + elif nT < 1.25*degreeT or nP < 1.25*degreeP: + logging.warning( + 'This Chebyshev fitting has few degrees of freedom and may not be '\ + 'accurate between data points. Consider increasing the number of '\ + 'temperature and pressure values in the fitting parameters.') # Set temperature and pressure ranges self.Tmin = (Tmin,"K") self.Tmax = (Tmax,"K") diff --git a/rmgpy/kinetics/chebyshevTest.py b/rmgpy/kinetics/chebyshevTest.py index 0b31c232a8..00447a9d87 100644 --- a/rmgpy/kinetics/chebyshevTest.py +++ b/rmgpy/kinetics/chebyshevTest.py @@ -36,7 +36,7 @@ import numpy from rmgpy.kinetics.chebyshev import Chebyshev - +from rmgpy.exceptions import KineticsError ################################################################################ class TestChebyshev(unittest.TestCase): @@ -153,7 +153,23 @@ def test_fitToData(self): for p in range(nP): kfit = chebyshev.getRateCoefficient(Tdata[t], Pdata[p]) * 1e6 self.assertAlmostEqual(kfit, kdata[t,p], delta=1e-4*kdata[t,p]) + + def test_fitToData2(self): + """ + Test the Chebyshev.fitToData() method throws error without enough degrees of freedom. + Here only 3 temperatures are given, but the polynomial desired has 6 parameters. + """ + Tdata = numpy.array([300,1200,2000]) + Pdata = numpy.array([1e5,3e5,1e6,3e7]) + nT = len(Tdata); nP = len(Pdata) + kdata = numpy.zeros((nT,nP)) + for t in range(nT): + for p in range(nP): + kdata[t,p] = self.chebyshev.getRateCoefficient(Tdata[t], Pdata[p]) + with self.assertRaises(KineticsError): + Chebyshev().fitToData(Tdata, Pdata, kdata, kunits="cm^3/(mol*s)", degreeT=12, degreeP=8, Tmin=300, Tmax=2000, Pmin=0.1, Pmax=10.) + def test_pickle(self): """ Test that a Chebyshev object can be pickled and unpickled with no loss diff --git a/rmgpy/molecule/atomtype.py b/rmgpy/molecule/atomtype.py index f88d5b9daf..02e5323ae2 100644 --- a/rmgpy/molecule/atomtype.py +++ b/rmgpy/molecule/atomtype.py @@ -68,8 +68,8 @@ class AtomType: `incrementLonePair` ``list`` The atom type(s) that result when the number of lone electron pairs is incremented `decrementLonePair` ``list`` The atom type(s) that result when the number of lone electron pairs is decremented - The following features are what are required in a given atomtype. Any int in the list is acceptable. - An empty list is a wildcard + The following features are what are required in a given atomtype. Any int in the list is acceptable. An empty list is a wildcard + ---------------------------------------------------------------------------- 'single' ''list'' The total number of single bonds on the atom 'allDouble' ''list'' The total number of double bonds on the atom 'rDouble' ''list'' The number of double bonds to any non-oxygen, nonsulfur @@ -235,25 +235,25 @@ def getFeatures(self): 'H', 'R!H', 'Val4','Val5','Val6','Val7', - 'He', + 'He','Ne','Ar', 'C','Ca','Cs','Csc','Cd','CO','CS','Cdd','Cdc','Ct','Cb','Cbf','C2s','C2sc','C2d','C2dc','C2tc', 'N','N0sc','N1s','N1sc','N1dc','N3s','N3sc','N3d','N3t','N3b','N5sc','N5dc','N5ddc','N5dddc','N5t','N5tc','N5b','N5bd', 'O','Oa','O0sc','O2s','O2sc','O2d','O4sc','O4dc','O4tc','O4b', - 'Ne', 'Si','Sis','Sid','Sidd','Sit','SiO','Sib','Sibf', 'S','Sa','S0sc','S2s','S2sc','S2d','S2dc','S2tc','S4s','S4sc','S4d','S4dd','S4dc','S4b','S4t','S4tdc','S6s','S6sc','S6d','S6dd','S6ddd','S6dc','S6t','S6td','S6tt','S6tdc', - 'Cl','Ar'] -) + 'Cl','Cl1s', + 'I','I1s']) + atomTypes['R!H'] = AtomType(label='R!H', generic=['R'], specific=[ - 'He', 'Val4','Val5','Val6','Val7', + 'He','Ne','Ar', 'C','Ca','Cs','Csc','Cd','CO','CS','Cdd','Cdc','Ct','Cb','Cbf','C2s','C2sc','C2d','C2dc','C2tc', 'N','N0sc','N1s','N1sc','N1dc','N3s','N3sc','N3d','N3t','N3b','N5sc','N5dc','N5ddc','N5dddc','N5t','N5tc','N5b','N5bd', 'O','Oa','O0sc','O2s','O2sc','O2d','O4sc','O4dc','O4tc','O4b', - 'Ne', 'Si','Sis','Sid','Sidd','Sit','SiO','Sib','Sibf', 'S','Sa','S0sc','S2s','S2sc','S2d','S2dc','S2tc','S4s','S4sc','S4d','S4dd','S4dc','S4b','S4t','S4tdc','S6s','S6sc','S6d','S6dd','S6ddd','S6dc','S6t','S6td','S6tt','S6tdc', - 'Cl','Ar']) + 'Cl','Cl1s', + 'I','I1s']) atomTypes['Val4'] = AtomType(label='Val4', generic=['R','R!H'], specific=[ 'C','Ca','Cs','Csc','Cd','CO','CS','Cdd','Cdc','Ct','Cb','Cbf','C2s','C2sc','C2d','C2dc','C2tc', @@ -267,11 +267,14 @@ def getFeatures(self): 'S','Sa','S0sc','S2s','S2sc','S2d','S2dc','S2tc','S4s','S4sc','S4d','S4dd','S4dc','S4b','S4t','S4tdc','S6s','S6sc','S6d','S6dd','S6ddd','S6dc','S6t','S6td','S6tt','S6tdc']) atomTypes['Val7'] = AtomType(label='Val7', generic=['R','R!H'], specific=[ - 'Cl']) + 'Cl','Cl1s', + 'I','I1s']) atomTypes['H' ] = AtomType('H', generic=['R'], specific=[]) -atomTypes['He' ] = AtomType('He', generic=['R','R!H'], specific=[]) +atomTypes['He' ] = AtomType('He', generic=['R','R!H'], specific=[]) +atomTypes['Ne' ] = AtomType('Ne', generic=['R','R!H'], specific=[]) +atomTypes['Ar' ] = AtomType('Ar', generic=['R','R!H'], specific=[]) atomTypes['C' ] = AtomType('C', generic=['R','R!H','Val4'], specific=['Ca','Cs','Csc','Cd','CO','CS','Cdd','Cdc','Ct','Cb','Cbf','C2s','C2sc','C2d','C2dc','C2tc'], single=[], allDouble=[], rDouble=[], oDouble=[], sDouble=[], triple=[], benzene=[], lonePairs=[], charge=[]) @@ -408,7 +411,6 @@ def getFeatures(self): single=[0], allDouble=[0], rDouble=[0], oDouble=[0], sDouble=[0], triple=[0], benzene=[2], lonePairs=[1], charge=[0]) # examples for S4b: Furane, Benzofurane, Benzo[c]thiophene, Oxazole... -atomTypes['Ne' ] = AtomType('Ne', generic=['R','R!H'], specific=[]) atomTypes['Si' ] = AtomType('Si', generic=['R','R!H','Val4'], specific=['Sis','Sid','Sidd','Sit','SiO','Sib','Sibf'], single=[], allDouble=[], rDouble=[], oDouble=[], sDouble=[], triple=[], benzene=[], lonePairs=[], charge=[]) atomTypes['Sis' ] = AtomType('Sis', generic=['R','R!H','Si','Val4'], specific=[], @@ -504,9 +506,15 @@ def getFeatures(self): single=[0,1,2,3,4], allDouble=[0,1,2], rDouble=[], oDouble=[], sDouble=[], triple=[1,2], benzene=[0], lonePairs=[0], charge=[-1,+1]) # *Composite atomType; examples for S6tdc: [SH2+]#[C-], [N-]=[S+]#N -atomTypes['Cl' ] = AtomType('Cl', generic=['R','R!H','Val7'], specific=[]) +atomTypes['Cl' ] = AtomType('Cl', generic=['R','R!H','Val7'], specific=['Cl1s']) +atomTypes['Cl1s'] = AtomType('Cl1s', generic=['R','R!H','Cl','Val7'], specific=[], + single=[0,1], allDouble=[0], rDouble=[], oDouble=[], sDouble=[], triple=[0], benzene=[0], lonePairs=[3], charge=[0]) +# examples for Cl1s: HCl, [Cl] -atomTypes['Ar' ] = AtomType('Ar', generic=['R','R!H'], specific=[]) +atomTypes['I' ] = AtomType('I', generic=['R','R!H','Val7'], specific=['I1s']) +atomTypes['I1s'] = AtomType('I1s', generic=['R','R!H','I','Val7'], specific=[], + single=[0,1], allDouble=[0], rDouble=[], oDouble=[], sDouble=[], triple=[0], benzene=[0], lonePairs=[3], charge=[0]) +# examples for I1s: HI, [I], IO, CH3I, I2 atomTypes['R' ].setActions(incrementBond=['R'], decrementBond=['R'], formBond=['R'], breakBond=['R'], incrementRadical=['R'], decrementRadical=['R'], incrementLonePair=['R'], decrementLonePair=['R']) atomTypes['R!H' ].setActions(incrementBond=['R!H'], decrementBond=['R!H'], formBond=['R!H'], breakBond=['R!H'], incrementRadical=['R!H'], decrementRadical=['R!H'], incrementLonePair=['R!H'], decrementLonePair=['R!H']) @@ -515,10 +523,11 @@ def getFeatures(self): atomTypes['Val6'].setActions(incrementBond=['Val6'], decrementBond=['Val6'], formBond=['Val6'], breakBond=['Val6'], incrementRadical=['Val6'], decrementRadical=['Val6'], incrementLonePair=['Val6'],decrementLonePair=['Val6']) atomTypes['Val7'].setActions(incrementBond=['Val7'], decrementBond=['Val7'], formBond=['Val7'], breakBond=['Val7'], incrementRadical=['Val7'], decrementRadical=['Val7'], incrementLonePair=['Val7'],decrementLonePair=['Val7']) - atomTypes['H' ].setActions(incrementBond=[], decrementBond=[], formBond=['H'], breakBond=['H'], incrementRadical=['H'], decrementRadical=['H'], incrementLonePair=[], decrementLonePair=[]) atomTypes['He' ].setActions(incrementBond=[], decrementBond=[], formBond=[], breakBond=[], incrementRadical=['He'], decrementRadical=['He'], incrementLonePair=[], decrementLonePair=[]) +atomTypes['Ne' ].setActions(incrementBond=[], decrementBond=[], formBond=[], breakBond=[], incrementRadical=['Ne'], decrementRadical=['Ne'], incrementLonePair=[], decrementLonePair=[]) +atomTypes['Ar' ].setActions(incrementBond=[], decrementBond=[], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[], incrementLonePair=[], decrementLonePair=[]) atomTypes['C' ].setActions(incrementBond=['C'], decrementBond=['C'], formBond=['C'], breakBond=['C'], incrementRadical=['C'], decrementRadical=['C'], incrementLonePair=['C'], decrementLonePair=['C']) atomTypes['Ca' ].setActions(incrementBond=[], decrementBond=[], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[], incrementLonePair=[], decrementLonePair=['C2s']) @@ -568,8 +577,6 @@ def getFeatures(self): atomTypes['O4tc'].setActions(incrementBond=[], decrementBond=['O4dc'], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[], incrementLonePair=['O2d'], decrementLonePair=[]) atomTypes['O4b' ].setActions(incrementBond=[], decrementBond=[], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[], incrementLonePair=[], decrementLonePair=[]) -atomTypes['Ne' ].setActions(incrementBond=[], decrementBond=[], formBond=[], breakBond=[], incrementRadical=['Ne'], decrementRadical=['Ne'], incrementLonePair=[], decrementLonePair=[]) - atomTypes['Si' ].setActions(incrementBond=['Si'], decrementBond=['Si'], formBond=['Si'], breakBond=['Si'], incrementRadical=['Si'], decrementRadical=['Si'], incrementLonePair=[], decrementLonePair=[]) atomTypes['Sis' ].setActions(incrementBond=['Sid','SiO'], decrementBond=[], formBond=['Sis'], breakBond=['Sis'], incrementRadical=['Sis'], decrementRadical=['Sis'], incrementLonePair=[], decrementLonePair=[]) atomTypes['Sid' ].setActions(incrementBond=['Sidd','Sit'], decrementBond=['Sis'], formBond=['Sid'], breakBond=['Sid'], incrementRadical=['Sid'], decrementRadical=['Sid'], incrementLonePair=[], decrementLonePair=[]) @@ -606,15 +613,16 @@ def getFeatures(self): atomTypes['S6tt'].setActions(incrementBond=[], decrementBond=['S6td','S6tdc'], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[], incrementLonePair=[], decrementLonePair=[]) atomTypes['S6tdc'].setActions(incrementBond=['S6td','S6tdc','S6tt'],decrementBond=['S6dc','S6tdc'],formBond=['S6tdc'],breakBond=['S6tdc'], incrementRadical=['S6tdc'],decrementRadical=['S6tdc'],incrementLonePair=['S4t','S4tdc'],decrementLonePair=[]) +atomTypes['Cl' ].setActions(incrementBond=[], decrementBond=[], formBond=['Cl'], breakBond=['Cl'], incrementRadical=['Cl'], decrementRadical=['Cl'], incrementLonePair=[], decrementLonePair=[]) +atomTypes['Cl1s'].setActions(incrementBond=[], decrementBond=[], formBond=['Cl1s'], breakBond=['Cl1s'], incrementRadical=['Cl1s'], decrementRadical=['Cl1s'], incrementLonePair=[], decrementLonePair=[]) -atomTypes['Cl' ].setActions(incrementBond=[], decrementBond=['Cl'], formBond=['Cl'], breakBond=['Cl'], incrementRadical=['Cl'], decrementRadical=['Cl'], incrementLonePair=[], decrementLonePair=[]) - -atomTypes['Ar' ].setActions(incrementBond=[], decrementBond=[], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[], incrementLonePair=[], decrementLonePair=[]) +atomTypes['I' ].setActions(incrementBond=[], decrementBond=[], formBond=['I'], breakBond=['I'], incrementRadical=['I'], decrementRadical=['I'], incrementLonePair=[], decrementLonePair=[]) +atomTypes['I1s'].setActions(incrementBond=[], decrementBond=[], formBond=['I1s'], breakBond=['I1s'], incrementRadical=['I1s'], decrementRadical=['I1s'], incrementLonePair=[], decrementLonePair=[]) #list of elements that do not have more specific atomTypes #these are ordered on priority of picking if we encounter a more general atomType for make allElements=['H', 'C', 'O', 'N', 'S', 'Si', 'Cl', 'Ne', 'Ar', 'He',] -nonSpecifics=['H', 'He', 'Ne', 'Cl', 'Ar',] +nonSpecifics=['H', 'He', 'Ne', 'Ar',] for atomType in atomTypes.values(): for items in [atomType.generic, atomType.specific, diff --git a/rmgpy/molecule/atomtypeTest.py b/rmgpy/molecule/atomtypeTest.py index 4cbd39fa75..b0c05fcfc6 100644 --- a/rmgpy/molecule/atomtypeTest.py +++ b/rmgpy/molecule/atomtypeTest.py @@ -486,6 +486,12 @@ def setUp(self): # 10 H u0 p0 {4,S} # 11 H u0 p0 {5,S}''') + self.mol73 = Molecule().fromAdjacencyList('''1 H u0 p0 c0 {2,S} + 2 Cl u0 p3 c0 {1,S}''') + + self.mol74 = Molecule().fromAdjacencyList('''1 H u0 p0 c0 {2,S} + 2 I u0 p3 c0 {1,S}''') + def atomType(self, mol, atomID): atom = mol.atoms[atomID] type = getAtomType(atom, mol.getBonds(atom)) @@ -604,7 +610,19 @@ def testSulfurTypes(self): self.assertEqual(self.atomType(self.mol36, 0), 'S6td') self.assertEqual(self.atomType(self.mol37, 1), 'S6tt') self.assertEqual(self.atomType(self.mol70, 0), 'S6tdc') - + + def testChlorineTypes(self): + """ + Test that getAtomType() returns appropriate chlorine atom types. + """ + self.assertEqual(self.atomType(self.mol73, 1), 'Cl1s') + + def testIodineTypes(self): + """ + Test that getAtomType() returns appropriate iodine atom types. + """ + self.assertEqual(self.atomType(self.mol74, 1), 'I1s') + def testOtherTypes(self): """ Test that getAtomType() returns appropriate types for other misc inerts. diff --git a/rmgpy/molecule/element.py b/rmgpy/molecule/element.py index 26c5ec9bc4..2e1318c148 100644 --- a/rmgpy/molecule/element.py +++ b/rmgpy/molecule/element.py @@ -52,17 +52,17 @@ class Element: """ A chemical element. The attributes are: - =========== =============== ================================================ - Attribute Type Description - =========== =============== ================================================ - `number` ``int`` The atomic number of the element - `symbol` ``str`` The symbol used for the element - `name` ``str`` The IUPAC name of the element - `mass` ``float`` The mass of the element in kg/mol - `covRadius` ``float`` Covalent bond radius in Angstrom - `isotope` ``int`` The isotope integer of the element - `chemkinName` ``str`` The chemkin compatible representation of the element - =========== =============== ================================================ + ============= =============== ================================================ + Attribute Type Description + ============= =============== ================================================ + `number` ``int`` The atomic number of the element + `symbol` ``str`` The symbol used for the element + `name` ``str`` The IUPAC name of the element + `mass` ``float`` The mass of the element in kg/mol + `covRadius` ``float`` Covalent bond radius in Angstrom + `isotope` ``int`` The isotope integer of the element + `chemkinName` ``str`` The chemkin compatible representation of the element + ============= =============== ================================================ This class is specifically for properties that all atoms of the same element share. Ideally there is only one instance of this class for each element. @@ -113,9 +113,9 @@ class PeriodicSystem(object): `lone_pairs`: the number of lone pairs an element has """ - valences = {'H': 1, 'He': 0, 'C': 4, 'N': 3, 'O': 2, 'Ne': 0, 'Si': 4, 'S': 2, 'Cl': 1, 'Ar': 0} - valence_electrons = {'H': 1, 'He': 2, 'C': 4, 'N': 5, 'O': 6, 'Ne': 8, 'Si': 4, 'S': 6, 'Cl': 7, 'Ar': 8} - lone_pairs = {'H': 0, 'He': 1, 'C': 0, 'N': 1, 'O': 2, 'Ne': 4, 'Si': 0, 'S': 2, 'Cl': 3, 'Ar': 4} + valences = {'H': 1, 'He': 0, 'C': 4, 'N': 3, 'O': 2, 'F': 1, 'Ne': 0, 'Si': 4, 'S': 2, 'Cl': 1, 'Ar': 0, 'I': 1} + valence_electrons = {'H': 1, 'He': 2, 'C': 4, 'N': 5, 'O': 6, 'F': 7, 'Ne': 8, 'Si': 4, 'S': 6, 'Cl': 7, 'Ar': 8, 'I': 7} + lone_pairs = {'H': 0, 'He': 1, 'C': 0, 'N': 1, 'O': 2, 'F': 3, 'Ne': 4, 'Si': 0, 'S': 2, 'Cl': 3, 'Ar': 4, 'I': 3} ################################################################################ diff --git a/rmgpy/molecule/generator.py b/rmgpy/molecule/generator.py index 124158369d..a10cae3a87 100644 --- a/rmgpy/molecule/generator.py +++ b/rmgpy/molecule/generator.py @@ -66,15 +66,16 @@ 'CH4O': 'CO', 'CO2': 'O=C=O', 'CO': '[C-]#[O+]', - 'C2H4': 'C=C', 'O2': 'O=O', 'C': '[C]', # for this to be in the "molecule" list it must be singlet with 2 lone pairs - 'SO2': 'O=S=O', - 'SO3': 'O=S(=O)=O', - 'H2SO4': 'OS(=O)(=O)O', + 'H2S': 'S', 'N2O': 'N#[N+][O-]', 'NH3': 'N', 'O3': '[O-][O+]=O', + 'Cl2': '[Cl][Cl]', + 'ClH': 'Cl', + 'I2': '[I][I]', + 'HI': 'I', } _known_smiles_radicals = { @@ -82,19 +83,23 @@ 'HO': '[OH]', 'C2H5': 'C[CH2]', 'O': '[O]', + 'S': '[S]', + 'N': '[N]', 'HO2': '[O]O', 'CH': '[CH]', + 'CH2': '[CH2]', 'H': '[H]', 'C': '[C]', # this, in the radical list, could be triplet or quintet. - #'CO2': it could be [O][C][O] or O=[C][O] - #'CO': '[C]=O', could also be [C][O] - #'C2H4': could be [CH3][CH] or [CH2][CH2] 'O2': '[O][O]', 'S2': '[S][S]', - 'SO': '[S][O]', - 'HSO3': 'OS(=O)[O]', + 'OS': '[S][O]', + 'HS': '[SH]', + 'H2N': '[NH2]', + 'HN': '[NH]', 'NO': '[N]=O', 'NO2': 'N(=O)[O]', + 'Cl': '[Cl]', + 'I': '[I]', } def toInChI(mol): diff --git a/rmgpy/molecule/group.pxd b/rmgpy/molecule/group.pxd index 7f2f0aedf7..88567ee57c 100644 --- a/rmgpy/molecule/group.pxd +++ b/rmgpy/molecule/group.pxd @@ -118,6 +118,9 @@ cdef class Group(Graph): cdef public short nitrogenCount cdef public short oxygenCount cdef public short sulfurCount + cdef public short chlorineCount + cdef public short iodineCount + cdef public short siliconCount cdef public short radicalCount cpdef addAtom(self, GroupAtom atom) diff --git a/rmgpy/molecule/group.py b/rmgpy/molecule/group.py index 0920205974..107765aa5d 100644 --- a/rmgpy/molecule/group.py +++ b/rmgpy/molecule/group.py @@ -627,6 +627,8 @@ def getOrderStr(self): values.append('T') elif value == 1.5: values.append('B') + elif value == 0: + values.append('H') else: raise TypeError('Bond order number {} is not hardcoded as a string'.format(value)) return values @@ -646,6 +648,8 @@ def setOrderStr(self, newOrder): values.append(3) elif value == 'B': values.append(1.5) + elif value == 'H': + values.append(0) else: # try to see if an float disguised as a string was input by mistake try: @@ -730,6 +734,21 @@ def isBenzene(self, wildcards = False): else: return abs(self.order[0]-1.5) <= 1e-9 and len(self.order) == 1 + def isHydrogenBond(self, wildcards = False): + """ + Return ``True`` if the bond represents a hydrogen bond or ``False`` if + not. If `wildcards` is ``False`` we return False anytime there is more + than one bond order, otherwise we return ``True`` if any of the options + are hydrogen bonds. + """ + if wildcards: + for order in self.order: + if abs(order) <= 1e-9: + return True + else: return False + else: + return abs(self.order[0]) <= 1e-9 and len(self.order) == 1 + def __changeBond(self, order): """ Update the bond group as a result of applying a CHANGE_BOND action, @@ -1098,18 +1117,26 @@ def updateFingerprint(self): isomorphism checks. """ cython.declare(atom=GroupAtom, atomType=AtomType) - cython.declare(carbon=AtomType, nitrogen=AtomType, oxygen=AtomType, sulfur=AtomType) - cython.declare(isCarbon=cython.bint, isNitrogen=cython.bint, isOxygen=cython.bint, isSulfur=cython.bint, radical=cython.int) - + cython.declare(carbon=AtomType, nitrogen=AtomType, oxygen=AtomType, sulfur=AtomType, chlorine=AtomType, + iodine=AtomType, silicon=AtomType) + cython.declare(isCarbon=cython.bint, isNitrogen=cython.bint, isOxygen=cython.bint, isSulfur=cython.bint, + isChlorine=cython.bint, isIodine=cython.bint, isSilicon=cython.bint, radical=cython.int) + carbon = atomTypes['C'] nitrogen = atomTypes['N'] oxygen = atomTypes['O'] sulfur = atomTypes['S'] + chlorine = atomTypes['Cl'] + iodine = atomTypes['I'] + silicon = atomTypes['Si'] self.carbonCount = 0 self.nitrogenCount = 0 self.oxygenCount = 0 self.sulfurCount = 0 + self.chlorineCount = 0 + self.iodineCount = 0 + self.siliconCount = 0 self.radicalCount = 0 for atom in self.vertices: if len(atom.atomType) == 1: @@ -1118,17 +1145,27 @@ def updateFingerprint(self): isNitrogen = atomType.equivalent(nitrogen) isOxygen = atomType.equivalent(oxygen) isSulfur = atomType.equivalent(sulfur) - if isCarbon and not isNitrogen and not isOxygen and not isSulfur: - self.carbonCount += 1 - elif isNitrogen and not isCarbon and not isOxygen and not isSulfur: - self.nitrogenCount += 1 - elif isOxygen and not isCarbon and not isNitrogen and not isSulfur: - self.oxygenCount += 1 - elif isSulfur and not isCarbon and not isNitrogen and not isOxygen: - self.sulfurCount += 1 - if len(atom.radicalElectrons) == 1: - radical = atom.radicalElectrons[0] - self.radicalCount += radical + isChlorine = atomType.equivalent(chlorine) + isIodine = atomType.equivalent(iodine) + isSilicon = atomType.equivalent(silicon) + sum_is_atom = isCarbon + isNitrogen + isOxygen + isSulfur + isChlorine + isIodine + isSilicon + if sum_is_atom == 1: + if isCarbon: + self.carbonCount += 1 + elif isNitrogen: + self.nitrogenCount += 1 + elif isOxygen: + self.oxygenCount += 1 + elif isSulfur: + self.sulfurCount += 1 + elif isChlorine: + self.chlorineCount += 1 + elif isIodine: + self.iodineCount += 1 + elif isSilicon: + self.siliconCount += 1 + if len(atom.radicalElectrons) >= 1: + self.radicalCount += atom.radicalElectrons[0] def isIsomorphic(self, other, initialMap=None): """ diff --git a/rmgpy/molecule/isomorphismTest.py b/rmgpy/molecule/isomorphismTest.py index c5a100050a..75dd788818 100644 --- a/rmgpy/molecule/isomorphismTest.py +++ b/rmgpy/molecule/isomorphismTest.py @@ -40,7 +40,7 @@ from rmgpy.molecule.molecule import Molecule from rmgpy.molecule.group import Group -molecule_atom_types = [ 'C', 'O', 'N', 'S', 'Si', 'Cl'] +molecule_atom_types = [ 'C', 'O', 'N', 'S', 'Si', 'Cl', 'I'] group_atomtypes = {} for item in create_atom_types() : @@ -140,6 +140,7 @@ def load_cases_molecule_atom_types(): ''' output = [] a_types = list(itertools.product(molecule_atom_types, repeat=2)) + uncharged_a_types = ['Cl','I'] unpaired_electrons = list(itertools.product(range(3), repeat=2)) cross_element_unpaired = list(itertools.product(a_types,unpaired_electrons)) for item in cross_element_unpaired: @@ -155,7 +156,10 @@ def load_cases_molecule_atom_types(): for now, only allow charges up to +1, not +2, +3, even if the unspecified valency allows for that. ''' - charges.append(range(min(val,1)+1)) + if el not in uncharged_a_types: + charges.append(range(min(val,1)+1)) + else: + charges.append((0,0)) charge_combos = list(itertools.product(charges[0],charges[1]))#cross product for both graphs for charge_combo in charge_combos:#combine charge tuple with the cross product of element and unpaired diff --git a/rmgpy/molecule/molecule.pxd b/rmgpy/molecule/molecule.pxd index 2b83565690..f9bf3a6506 100644 --- a/rmgpy/molecule/molecule.pxd +++ b/rmgpy/molecule/molecule.pxd @@ -61,7 +61,15 @@ cdef class Atom(Vertex): cpdef bint isOxygen(self) + cpdef bint isFluorine(self) + + cpdef bint isSilicon(self) + cpdef bint isSulfur(self) + + cpdef bint isChlorine(self) + + cpdef bint isIodine(self) cpdef incrementRadical(self) diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index 115e562875..4c722abc20 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -43,6 +43,7 @@ import urllib from collections import OrderedDict import itertools +from copy import deepcopy import element as elements try: @@ -51,6 +52,7 @@ pass from .graph import Vertex, Edge, Graph, getVertexConnectivityValue import rmgpy.molecule.group as gr +from rmgpy.molecule.pathfinder import find_shortest_path from .atomtype import AtomType, atomTypes, getAtomType, AtomTypeError import rmgpy.constants as constants import rmgpy.molecule.parser as parser @@ -305,20 +307,41 @@ def isOxygen(self): """ return self.element.number == 8 + def isFluorine(self): + """ + Return ``True`` if the atom represents a fluorine atom or ``False`` if + not. + """ + return self.element.number == 9 + def isSilicon(self): """ - Return ``True`` if the atom represents an silicon atom or ``False`` if + Return ``True`` if the atom represents a silicon atom or ``False`` if not. """ return self.element.number == 14 def isSulfur(self): """ - Return ``True`` if the atom represents an sulfur atom or ``False`` if + Return ``True`` if the atom represents a sulfur atom or ``False`` if not. """ return self.element.number == 16 + def isChlorine(self): + """ + Return ``True`` if the atom represents a chlorine atom or ``False`` if + not. + """ + return self.element.number == 17 + + def isIodine(self): + """ + Return ``True`` if the atom represents an iodine atom or ``False`` if + not. + """ + return self.element.number == 53 + def incrementRadical(self): """ Update the atom pattern as a result of applying a GAIN_RADICAL action, @@ -512,6 +535,8 @@ def getOrderStr(self): return 'D' elif self.isTriple(): return 'T' + elif self.isHydrogenBond(): + return 'H' else: raise ValueError("Bond order {} does not have string representation.".format(self.order)) @@ -527,6 +552,8 @@ def setOrderStr(self, newOrder): self.order = 3 elif newOrder == 'B': self.order = 1.5 + elif newOrder == 'H': + self.order = 0 else: # try to see if an float disguised as a string was input by mistake try: @@ -564,7 +591,7 @@ def copy(self): def isOrder(self, otherOrder): """ - Return ``True`` if the bond represents a single bond or ``False`` if + Return ``True`` if the bond is of order otherOrder or ``False`` if not. This compares floats that takes into account floating point error NOTE: we can replace the absolute value relation with math.isclose when @@ -600,7 +627,14 @@ def isBenzene(self): not. """ return self.isOrder(1.5) - + + def isHydrogenBond(self): + """ + Return ``True`` if the bond represents a hydrogen bond or ``False`` if + not. + """ + return self.isOrder(0) + def incrementOrder(self): """ Update the bond as a result of applying a CHANGE_BOND action to @@ -1455,7 +1489,81 @@ def toAdjacencyList(self, label='', removeH=False, removeLonePairs=False, oldSty from .adjlist import toAdjacencyList result = toAdjacencyList(self.vertices, self.multiplicity, label=label, group=False, removeH=removeH, removeLonePairs=removeLonePairs, oldStyle=oldStyle) return result + + def find_H_bonds(self): + """ + generates a list of (new-existing H bonds ignored) possible Hbond coordinates [(i1,j1),(i2,j2),...] where i and j values + correspond to the indexes of the atoms involved, Hbonds are allowed if they meet + the following constraints: + 1) between a H and [O,N] atoms + 2) the hydrogen is covalently bonded to an O or N + 3) the Hydrogen bond must complete a ring with at least 5 members + 4) An atom can only be hydrogen bonded to one other atom + """ + potBonds = [] + + ONatoms = [a for a in self.atoms if a.isOxygen() or a.isNitrogen()] + ONinds = [n for n,a in enumerate(self.atoms) if a.isOxygen() or a.isNitrogen()] + + for i,atm1 in enumerate(self.atoms): + if atm1.atomType.label == 'H': + atm_covs = [q for q in atm1.bonds.keys()] + if len(atm_covs) > 1: #H is already H bonded + continue + else: + atm_cov = atm_covs[0] + if (atm_cov.isOxygen() or atm_cov.isNitrogen()): #this H can be H-bonded + for k,atm2 in enumerate(ONatoms): + if all([q.order != 0 for q in atm2.bonds.values()]): #atm2 not already H bonded + dist = len(find_shortest_path(atm1,atm2))-1 + if dist > 3: + j = ONinds[k] + potBonds.append((i,j)) + return potBonds + + def generate_H_bonded_structures(self): + """ + generates a list of Hbonded molecular structures in addition to the + constraints on Hydrogen bonds applied in the find_H_Bonds function + the generated structures are constrained to: + + 1) An atom can only be hydrogen bonded to one other atom + 2) Only two H-bonds can exist in a given molecule + + the second is done to avoid explosive growth in the number of + structures as without this constraint the number of possible + structures grows 2^n where n is the number of possible H-bonds + """ + structs = [] + Hbonds = self.find_H_bonds() + for i,bd1 in enumerate(Hbonds): + molc = deepcopy(self) + molc.addBond(Bond(molc.atoms[bd1[0]],molc.atoms[bd1[1]],order=0)) + structs.append(molc) + for j,bd2 in enumerate(Hbonds): + if j 0 and densStates[r,s] > 0: k[r,s] = sumStates[r,s] / densStates[r,s] * dE - + logging.debug('Finished applying RRKM for path transition state {}'.format(transitionState)) + logging.debug('The rate constant is found to be {}'.format(k)) + return k @cython.boundscheck(False) @@ -305,7 +309,7 @@ def applyInverseLaplaceTransformMethod(transitionState, phi0 = numpy.zeros(Ngrains, numpy.float64) for r in range(Ngrains): E = Elist[r] - Elist[0] - Ea - if E > 0: + if E > 1: phi0[r] = (E/R)**(n-1.0) phi0 = phi0 * (dE / R) / scipy.special.gamma(n) # Evaluate the convolution @@ -318,6 +322,8 @@ def applyInverseLaplaceTransformMethod(transitionState, else: raise Exception('Unable to use inverse Laplace transform method for non-Arrhenius kinetics or for n < 0.') + logging.debug('Finished applying inverse lapace transform for path transition state {}'.format(transitionState)) + logging.debug('The rate constant is found to be {}'.format(k)) return k @@ -383,6 +389,7 @@ def fitInterpolationModel(reaction, Tlist, Plist, K, model, Tmin, Tmax, Pmin, Pm logRMS = sqrt(logRMS / len(Tlist) / len(Plist)) if logRMS > 0.5: logging.warning('RMS error for k(T,P) fit = {0:g} for reaction {1}.'.format(logRMS, reaction)) - + logging.debug('Finished fitting model for path reaction {}'.format(reaction)) + logging.debug('The kinetics fit is {0!r}'.format(kinetics)) return kinetics diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py index 30bc1ffea8..8aa6863283 100644 --- a/rmgpy/reaction.py +++ b/rmgpy/reaction.py @@ -366,8 +366,8 @@ def matchesSpecies(self, reactants, products=None): and products of this reaction. Both directions are checked. Args: - reactants List of Species required on one side of the reaction - products List of Species required on the other side (optional) + reactants (list): Species required on one side of the reaction + products (list, optional): Species required on the other side """ # Check forward direction if _isomorphicSpeciesList(self.reactants, reactants): @@ -916,24 +916,30 @@ def generatePairs(self): productOxygens = [sum([1 for atom in product.molecule[0].atoms if atom.isOxygen()]) for product in products ] reactantNitrogens = [sum([1 for atom in reactant.molecule[0].atoms if atom.isNitrogen()]) for reactant in reactants] productNitrogens = [sum([1 for atom in product.molecule[0].atoms if atom.isNitrogen()]) for product in products ] + reactantSilicons = [sum([1 for atom in reactant.molecule[0].atoms if atom.isSilicon()]) for reactant in reactants] + productSilicons = [sum([1 for atom in product.molecule[0].atoms if atom.isSilicon()]) for product in products ] reactantSulfurs = [sum([1 for atom in reactant.molecule[0].atoms if atom.isSulfur()]) for reactant in reactants] productSulfurs = [sum([1 for atom in product.molecule[0].atoms if atom.isSulfur()]) for product in products ] + reactantChlorines = [sum([1 for atom in reactant.molecule[0].atoms if atom.isChlorine()]) for reactant in reactants] + productChlorines = [sum([1 for atom in product.molecule[0].atoms if atom.isChlorine()]) for product in products ] + reactantIodines = [sum([1 for atom in reactant.molecule[0].atoms if atom.isChlorine()]) for reactant in reactants] + productIodines = [sum([1 for atom in product.molecule[0].atoms if atom.isChlorine()]) for product in products ] # Sort the reactants and products by C/O/N/S numbers - reactants = [(carbon, oxygen, nitrogen, sulfur, reactant) for carbon, oxygen, nitrogen, sulfur, reactant - in zip(reactantCarbons,reactantOxygens,reactantNitrogens,reactantSulfurs,reactants)] + reactants = [(carbon, oxygen, nitrogen, silicon, sulfur, chlorine, iodine, reactant) for carbon, oxygen, nitrogen, silicon, sulfur, chlorine, iodine, reactant + in zip(reactantCarbons,reactantOxygens,reactantNitrogens,reactantSilicons,reactantSulfurs,reactantChlorines, reactantIodines, reactants)] reactants.sort() - products = [(carbon, oxygen, nitrogen, sulfur, product) for carbon, oxygen, nitrogen, sulfur, product - in zip(productCarbons,productOxygens,productNitrogens,productSulfurs,products)] + products = [(carbon, oxygen, nitrogen, silicon, sulfur, chlorine, iodine, product) for carbon, oxygen, nitrogen, silicon, sulfur, chlorine, iodine, product + in zip(productCarbons,productOxygens,productNitrogens,productSilicons,productSulfurs,productChlorines, productIodines, products)] products.sort() while len(reactants) > 1 and len(products) > 1: - self.pairs.append((reactants[-1][4], products[-1][4])) + self.pairs.append((reactants[-1][-1], products[-1][-1])) reactants.pop() products.pop() for reactant in reactants: for product in products: - self.pairs.append((reactant[4], product[4])) + self.pairs.append((reactant[-1], product[-1])) def draw(self, path): """ @@ -1071,11 +1077,11 @@ def ensure_species(self, reactant_resonance=False, product_resonance=True): return None # obtain species with all resonance isomers if self.isForward: - self.reactants = ensure_species(self.reactants, resonance=reactant_resonance, keepIsomorphic=True) - self.products = ensure_species(self.products, resonance=product_resonance, keepIsomorphic=True) + ensure_species(self.reactants, resonance=reactant_resonance, keepIsomorphic=True) + ensure_species(self.products, resonance=product_resonance, keepIsomorphic=True) else: - self.reactants = ensure_species(self.reactants, resonance=product_resonance, keepIsomorphic=True) - self.products = ensure_species(self.products, resonance=reactant_resonance, keepIsomorphic=True) + ensure_species(self.reactants, resonance=product_resonance, keepIsomorphic=True) + ensure_species(self.products, resonance=reactant_resonance, keepIsomorphic=True) # convert reaction.pairs object to species if self.pairs: diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py index 57ace7cba8..ac9fa1587b 100644 --- a/rmgpy/rmg/input.py +++ b/rmgpy/rmg/input.py @@ -104,6 +104,10 @@ def database( def species(label, structure, reactive=True): logging.debug('Found {0} species "{1}" ({2})'.format('reactive' if reactive else 'nonreactive', label, structure.toSMILES())) + + if '+' in label: + raise InputError('species {0} label cannot include a + sign'.format(label)) + spec, isNew = rmg.reactionModel.makeNewSpecies(structure, label=label, reactive=reactive) if not isNew: raise InputError("Species {0} is a duplicate of {1}. Species in input file must be unique".format(label,spec.label)) diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index f23b796164..0366ce37a5 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -70,8 +70,9 @@ from rmgpy.qm.main import QMDatabaseWriter from rmgpy.stats import ExecutionStatsWriter from rmgpy.thermo.thermoengine import submit -from rmgpy.tools.sensitivity import plotSensitivity +from rmgpy.tools.simulate import plot_sensitivity from cantera import ck2cti +from rmgpy.exceptions import CoreError ################################################################################ solvent = None @@ -785,7 +786,7 @@ def execute(self, **kwargs): simulatorSettings = self.simulatorSettingsList[-1], ) - plotSensitivity(self.outputDirectory, index, reactionSystem.sensitiveSpecies) + plot_sensitivity(self.outputDirectory, index, reactionSystem.sensitiveSpecies) # generate Cantera files chem.cti & chem_annotated.cti in a designated `cantera` output folder try: @@ -793,7 +794,9 @@ def execute(self, **kwargs): self.generateCanteraFiles(os.path.join(self.outputDirectory, 'chemkin', 'chem_annotated.inp')) except EnvironmentError: logging.error('Could not generate Cantera files due to EnvironmentError. Check read\write privileges in output directory.') - + + self.check_model() + # Write output file logging.info('') logging.info('MODEL GENERATION COMPLETED') @@ -804,6 +807,34 @@ def execute(self, **kwargs): self.finish() + def check_model(self): + """ + Run checks on the RMG model + """ + logging.info('Performing final model checks...') + + # Check that no two species in core or edge are isomorphic + for i, spc in enumerate(self.reactionModel.core.species): + for j in xrange(i): + spc2 = self.reactionModel.core.species[j] + if spc.isIsomorphic(spc2): + raise CoreError( + 'Although the model has completed, species {0} is isomorphic to species {1} in the core. ' + 'Please open an issue on GitHub with the following output:' + '\n{2}\n{3}'.format(spc.label, spc2.label, spc.toAdjacencyList(), spc2.toAdjacencyList()) + ) + + for i, spc in enumerate(self.reactionModel.edge.species): + for j in xrange(i): + spc2 = self.reactionModel.edge.species[j] + if spc.isIsomorphic(spc2): + logging.warning( + 'Species {0} is isomorphic to species {1} in the edge. This does not affect ' + 'the generated model. If you would like to report this to help make RMG better ' + 'please open a GitHub issue with the following output:' + '\n{2}\n{3}'.format(spc.label, spc2.label, spc.toAdjacencyList(), spc2.toAdjacencyList()) + ) + def makeSeedMech(self,firstTime=False): """ causes RMG to make a seed mechanism out of the current chem_annotated.inp and species_dictionary.txt @@ -845,31 +876,7 @@ def makeSeedMech(self,firstTime=False): oldLabels = self.makeSpeciesLabelsIndependent(speciesList) edgeOldLabels = self.makeSpeciesLabelsIndependent(edgeSpeciesList) - # load thermo library entries - thermoLibrary = ThermoLibrary(name=name) - for i,species in enumerate(speciesList): - if species.thermo: - thermoLibrary.loadEntry(index = i + 1, - label = species.label, - molecule = species.molecule[0].toAdjacencyList(), - thermo = species.thermo, - shortDesc = species.thermo.comment - ) - else: - logging.warning('Species {0} did not contain any thermo data and was omitted from the thermo library.'.format(str(species))) - - edgeThermoLibrary = ThermoLibrary(name=name+'_edge') - for i,species in enumerate(edgeSpeciesList): - if species.thermo: - edgeThermoLibrary.loadEntry(index = i + 1, - label = species.label, - molecule = species.molecule[0].toAdjacencyList(), - thermo = species.thermo, - shortDesc = species.thermo.comment - ) - else: - logging.warning('Species {0} did not contain any thermo data and was omitted from the edge thermo library.'.format(str(species))) - + # load kinetics library entries kineticsLibrary = KineticsLibrary(name=name,autoGenerated=True) kineticsLibrary.entries = {} @@ -881,17 +888,16 @@ def makeSeedMech(self,firstTime=False): item = reaction, data = reaction.kinetics, ) - try: + + if 'rate rule' in reaction.kinetics.comment: + entry.longDesc = reaction.kinetics.comment + elif hasattr(reaction,'library') and reaction.library: entry.longDesc = 'Originally from reaction library: ' + reaction.library + "\n" + reaction.kinetics.comment - except AttributeError: + else: entry.longDesc = reaction.kinetics.comment + kineticsLibrary.entries[i+1] = entry - # Mark as duplicates where there are mixed pressure dependent and non-pressure dependent duplicate kinetics - # Even though CHEMKIN does not require a duplicate flag, RMG needs it. - # Using flag markDuplicates = True - - # load kinetics library entries edgeKineticsLibrary = KineticsLibrary(name=name+'_edge',autoGenerated=True) edgeKineticsLibrary.entries = {} @@ -903,16 +909,11 @@ def makeSeedMech(self,firstTime=False): data = reaction.kinetics, ) try: - entry.longDesc = 'Originally from reaction library: ' + reaction.library + "\n" + reaction.kinetics.comment - except AttributeError: - entry.longDesc = reaction.kinetics.comment + entry.longDesc = 'Originally from reaction library: ' + reaction.library + "\n" + reaction.kinetics.comment + except AttributeError: + entry.longDesc = reaction.kinetics.comment edgeKineticsLibrary.entries[i+1] = entry - # Mark as duplicates where there are mixed pressure dependent and non-pressure dependent duplicate kinetics - # Even though CHEMKIN does not require a duplicate flag, RMG needs it. - # Using flag markDuplicates = True - - #save in database if self.saveSeedToDatabase: databaseDirectory = settings['database.directory'] @@ -920,7 +921,6 @@ def makeSeedMech(self,firstTime=False): os.makedirs(os.path.join(databaseDirectory, 'kinetics', 'libraries',name)) except: pass - thermoLibrary.save(os.path.join(databaseDirectory, 'thermo' ,'libraries', name + '.py')) kineticsLibrary.save(os.path.join(databaseDirectory, 'kinetics', 'libraries', name, 'reactions.py')) kineticsLibrary.saveDictionary(os.path.join(databaseDirectory, 'kinetics', 'libraries', name, 'dictionary.txt')) @@ -928,16 +928,13 @@ def makeSeedMech(self,firstTime=False): os.makedirs(os.path.join(databaseDirectory, 'kinetics', 'libraries',name+'_edge')) except: pass - edgeThermoLibrary.save(os.path.join(databaseDirectory, 'thermo' ,'libraries', name +'_edge'+'.py')) edgeKineticsLibrary.save(os.path.join(databaseDirectory, 'kinetics', 'libraries', name+'_edge', 'reactions.py')) edgeKineticsLibrary.saveDictionary(os.path.join(databaseDirectory, 'kinetics', 'libraries', name+'_edge', 'dictionary.txt')) #save in output directory - thermoLibrary.save(os.path.join(seedDir, name + '.py')) kineticsLibrary.save(os.path.join(seedDir, name, 'reactions.py')) kineticsLibrary.saveDictionary(os.path.join(seedDir, name, 'dictionary.txt')) - edgeThermoLibrary.save(os.path.join(seedDir, name + '_edge'+ '.py')) edgeKineticsLibrary.save(os.path.join(seedDir, name+'_edge', 'reactions.py')) edgeKineticsLibrary.saveDictionary(os.path.join(seedDir, name+'_edge', 'dictionary.txt')) @@ -959,12 +956,19 @@ def makeSpeciesLabelsIndependent(self, species): for spec in species: oldLabels.append(spec.label) duplicate_index = 1 - potential_label = spec.label + if '+' in spec.label: + L = spec.molecule[0].getFormula() + else: + L = spec.label + potential_label = L while potential_label in labels: duplicate_index += 1 - potential_label = spec.label + '-{}'.format(duplicate_index) + potential_label = L + '-{}'.format(duplicate_index) + spec.label = potential_label labels.add(potential_label) + + return oldLabels ################################################################################ diff --git a/rmgpy/rmg/mainTest.py b/rmgpy/rmg/mainTest.py index b5b4980683..034dfce771 100644 --- a/rmgpy/rmg/mainTest.py +++ b/rmgpy/rmg/mainTest.py @@ -50,9 +50,7 @@ def setUpClass(cls): cls.outputDir = 'output' cls.databaseDirectory = settings['database.directory'] - cls.seedThermo = os.path.join(cls.databaseDirectory, 'thermo', 'libraries', 'testSeed.py') cls.seedKinetics = os.path.join(cls.databaseDirectory, 'kinetics', 'libraries', 'testSeed') - cls.seedThermoEdge = os.path.join(cls.databaseDirectory, 'thermo', 'libraries', 'testSeed_edge.py') cls.seedKineticsEdge = os.path.join(cls.databaseDirectory, 'kinetics', 'libraries', 'testSeed_edge') os.mkdir(os.path.join(cls.testDir, cls.outputDir)) @@ -73,9 +71,7 @@ def tearDownClass(cls): shutil.rmtree(os.path.join(cls.testDir, cls.outputDir)) # Delete the seed libraries created in database - os.remove(cls.seedThermo) shutil.rmtree(cls.seedKinetics) - os.remove(cls.seedThermoEdge) shutil.rmtree(cls.seedKineticsEdge) def testRMGExecute(self): @@ -95,7 +91,6 @@ def testRMGSeedMechanismCreation(self): seedDir = os.path.join(self.testDir, self.outputDir, 'seed') self.assertTrue(os.path.exists) - self.assertTrue(os.path.exists(os.path.join(seedDir, self.rmg.name + '.py'))) # thermo library made self.assertTrue(os.path.exists(os.path.join(seedDir, self.rmg.name))) # kinetics library folder made self.assertTrue(os.path.exists(os.path.join(seedDir, self.rmg.name, 'dictionary.txt'))) # dictionary file made @@ -108,7 +103,6 @@ def testRMGSeedEdgeMechanismCreation(self): name = self.rmg.name + '_edge' - self.assertTrue(os.path.exists(os.path.join(seedDir, name + '.py'))) # thermo library made self.assertTrue(os.path.exists(os.path.join(seedDir, name))) # kinetics library folder made self.assertTrue(os.path.exists(os.path.join(seedDir, name, 'dictionary.txt'))) # dictionary file made @@ -116,12 +110,10 @@ def testRMGSeedEdgeMechanismCreation(self): def testRMGSeedLibraryCreation(self): """Test that seed mechanisms are created in the correct database locations.""" - self.assertTrue(os.path.exists(self.seedThermo)) self.assertTrue(os.path.exists(self.seedKinetics)) def testRMGSeedEdgeLibraryCreation(self): """Test that edge seed mechanisms are created in the correct database locations.""" - self.assertTrue(os.path.exists(self.seedThermo)) self.assertTrue(os.path.exists(self.seedKinetics)) def testRMGSeedWorks(self): @@ -132,7 +124,7 @@ def testRMGSeedWorks(self): # Load the seed libraries into the database self.rmg.database.load( path=self.databaseDirectory, - thermoLibraries=['testSeed', 'testSeed_edge'], + thermoLibraries=[], reactionLibraries=['testSeed', 'testSeed_edge'], seedMechanisms=['testSeed', 'testSeed_edge'], kineticsFamilies='default', diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index 1b1c769506..0c3b6e70ea 100644 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -1401,6 +1401,8 @@ def addSeedMechanismToCore(self, seedMechanism, react=False): database = rmgpy.data.rmg.database libraryNames = database.kinetics.libraries.keys() + familyNames = database.kinetics.families.keys() + path = os.path.join(settings['database.directory'],'kinetics','libraries') from rmgpy.rmg.input import rmg @@ -1414,11 +1416,18 @@ def addSeedMechanismToCore(self, seedMechanism, react=False): seedMechanism = database.kinetics.libraries[seedMechanism] rxns = seedMechanism.getLibraryReactions() + for rxn in rxns: if isinstance(rxn,LibraryReaction) and not (rxn.library in libraryNames): #if one of the reactions in the library is from another library load that library database.kinetics.libraryOrder.append((rxn.library,'Internal')) database.kinetics.loadLibraries(path=path,libraries=[rxn.library]) libraryNames = database.kinetics.libraries.keys() + if isinstance(rxn,TemplateReaction) and not (rxn.family in familyNames): + logging.warning('loading reaction {0} originally from family {1} as a library reaction'.format(str(rxn),rxn.family)) + rxn = LibraryReaction(reactants=rxn.reactants[:], products=rxn.products[:], + library=seedMechanism.name, specificCollider=rxn.specificCollider, kinetics=rxn.kinetics, duplicate=rxn.duplicate, + reversible=rxn.reversible + ) r, isNew = self.makeNewReaction(rxn) # updates self.newSpeciesList and self.newReactionlist if not isNew: logging.info("This library reaction was not new: {0}".format(rxn)) @@ -1473,6 +1482,7 @@ def addReactionLibraryToEdge(self, reactionLibrary): database = rmgpy.data.rmg.database libraryNames = database.kinetics.libraries.keys() + familyNames = database.kinetics.families.keys() path = os.path.join(settings['database.directory'],'kinetics','libraries') from rmgpy.rmg.input import rmg @@ -1492,6 +1502,12 @@ def addReactionLibraryToEdge(self, reactionLibrary): database.kinetics.libraryOrder.append((rxn.library,'Internal')) database.kinetics.loadLibraries(path=path,libraries=[rxn.library]) libraryNames = database.kinetics.libraries.keys() + if isinstance(rxn,TemplateReaction) and not (rxn.family in familyNames): + logging.warning('loading reaction {0} originally from family {1} as a library reaction'.format(str(rxn),rxn.family)) + rxn = LibraryReaction(reactants=rxn.reactants[:], products=rxn.products[:], + library=reactionLibrary.name, specificCollider=rxn.specificCollider, kinetics=rxn.kinetics, duplicate=rxn.duplicate, + reversible=rxn.reversible + ) r, isNew = self.makeNewReaction(rxn) # updates self.newSpeciesList and self.newReactionlist if not isNew: logging.info("This library reaction was not new: {0}".format(rxn)) diff --git a/rmgpy/solver/__init__.py b/rmgpy/solver/__init__.py index ab7a795888..b2b0cd5fae 100644 --- a/rmgpy/solver/__init__.py +++ b/rmgpy/solver/__init__.py @@ -30,3 +30,4 @@ from .base import ReactionSystem, TerminationTime, TerminationConversion from .simple import SimpleReactor +from .liquid import LiquidReactor diff --git a/rmgpy/solver/baseTest.py b/rmgpy/solver/baseTest.py index bc2d025ce5..5268893084 100644 --- a/rmgpy/solver/baseTest.py +++ b/rmgpy/solver/baseTest.py @@ -55,7 +55,7 @@ def setUp(self): chemkinFile = os.path.join(folder, 'chemkin/chem.inp') spc_dict = os.path.join(folder, 'chemkin/species_dictionary.txt') - self.rmg = loadRMGPyJob(inputFile, chemkinFile, spc_dict, generateImages=False) + self.rmg = loadRMGPyJob(inputFile, chemkinFile, spc_dict, generateImages=False, checkDuplicates=False) def testSurfaceInitialization(self): """ diff --git a/rmgpy/solver/simple.pyx b/rmgpy/solver/simple.pyx index ddda83669f..157e162105 100644 --- a/rmgpy/solver/simple.pyx +++ b/rmgpy/solver/simple.pyx @@ -178,13 +178,18 @@ cdef class SimpleReactor(ReactionSystem): def calculate_effective_pressure(self, rxn): """ Computes the effective pressure for a reaction as: - Peff = P * sum(yi * effi / sum(y)) + + .. math:: P_{eff} = P * \\sum_i \\frac{y_i * eff_i}{\\sum_j y_j} + with: - P the pressure of the reactor, - y the array of initial moles of the core species + or as: - Peff = P * y_specificCollider / sum(y) - if a secificCollider is mentioned + + .. math:: P_{eff} = \\frac{P * y_{specificCollider}}{\\sum_j y_j} + + if a specificCollider is mentioned. """ y0_coreSpecies = self.y0[:self.numCoreSpecies] diff --git a/rmgpy/species.py b/rmgpy/species.py index 701f0fea0d..0ac9f0d9d8 100644 --- a/rmgpy/species.py +++ b/rmgpy/species.py @@ -50,6 +50,7 @@ import rmgpy.quantity as quantity +from rmgpy.molecule.molecule import Atom, Bond, Molecule from rmgpy.pdep import SingleExponentialDown from rmgpy.statmech.conformer import Conformer from rmgpy.thermo import Wilhoit, NASA, ThermoData diff --git a/rmgpy/statmech/conformer.pyx b/rmgpy/statmech/conformer.pyx index 43f700441f..a436b2e0fc 100644 --- a/rmgpy/statmech/conformer.pyx +++ b/rmgpy/statmech/conformer.pyx @@ -337,14 +337,18 @@ cdef class Conformer: Finally, the reduced moment of inertia is evaluated from the moment of inertia of each top via the formula (I1*I2)/(I1+I2). - option corresponds to 3 possible ways of calculating the internal reduced moment of inertia + Option corresponds to 3 possible ways of calculating the internal reduced moment of inertia as discussed in East and Radom [2] - option = 1 -> moments of inertia of each rotating group calculated about the axis containing the twisting bond - option = 2 (unimplemented) -> each moment of inertia of each rotating group is calculated about an axis parallel to the - twisting bond and passing through its center of mass - option = 3 -> moments of inertia of each rotating group calculated about the axis passing through the - centers of mass of both groups + +----------+---------------------------------------------------------------------------------------------------+ + |option = 1|moments of inertia of each rotating group calculated about the axis containing the twisting bond | + +----------+---------------------------------------------------------------------------------------------------+ + |option = 2|(unimplemented) each moment of inertia of each rotating group is calculated about an axis parallel | + | |to the twisting bond and passing through its center of mass | + +----------+---------------------------------------------------------------------------------------------------+ + |option = 3|moments of inertia of each rotating group calculated about the axis passing through the | + | |centers of mass of both groups | + +----------+---------------------------------------------------------------------------------------------------+ .. math:: \\frac{1}{I^{(2,option)}} = \\frac{1}{I_1} + \\frac{1}{I_2} diff --git a/rmgpy/tools/data/sens/liquid/chem.inp b/rmgpy/tools/data/sim/liquid/chem.inp similarity index 100% rename from rmgpy/tools/data/sens/liquid/chem.inp rename to rmgpy/tools/data/sim/liquid/chem.inp diff --git a/rmgpy/tools/data/sens/liquid/input.py b/rmgpy/tools/data/sim/liquid/input.py similarity index 100% rename from rmgpy/tools/data/sens/liquid/input.py rename to rmgpy/tools/data/sim/liquid/input.py diff --git a/rmgpy/tools/data/sens/liquid/species_dictionary.txt b/rmgpy/tools/data/sim/liquid/species_dictionary.txt similarity index 100% rename from rmgpy/tools/data/sens/liquid/species_dictionary.txt rename to rmgpy/tools/data/sim/liquid/species_dictionary.txt diff --git a/rmgpy/tools/data/sens/simple/chem.inp b/rmgpy/tools/data/sim/simple/chem.inp similarity index 100% rename from rmgpy/tools/data/sens/simple/chem.inp rename to rmgpy/tools/data/sim/simple/chem.inp diff --git a/rmgpy/tools/data/sens/simple/input.py b/rmgpy/tools/data/sim/simple/input.py similarity index 100% rename from rmgpy/tools/data/sens/simple/input.py rename to rmgpy/tools/data/sim/simple/input.py diff --git a/rmgpy/tools/data/sens/simple/species_dictionary.txt b/rmgpy/tools/data/sim/simple/species_dictionary.txt similarity index 100% rename from rmgpy/tools/data/sens/simple/species_dictionary.txt rename to rmgpy/tools/data/sim/simple/species_dictionary.txt diff --git a/rmgpy/tools/fluxdiagram.py b/rmgpy/tools/fluxdiagram.py index d940543044..875a574a17 100644 --- a/rmgpy/tools/fluxdiagram.py +++ b/rmgpy/tools/fluxdiagram.py @@ -442,7 +442,7 @@ def loadChemkinOutput(outputFile, reactionModel): ################################################################################ def createFluxDiagram(inputFile, chemkinFile, speciesDict, savePath=None, speciesPath=None, java=False, settings=None, - chemkinOutput='', centralSpecies=None, diffusionLimited=True): + chemkinOutput='', centralSpecies=None, diffusionLimited=True, checkDuplicates=True): """ Generates the flux diagram based on a condition 'inputFile', chemkin.inp chemkinFile, a speciesDict txt file, plus an optional chemkinOutput file. @@ -454,7 +454,8 @@ def createFluxDiagram(inputFile, chemkinFile, speciesDict, savePath=None, specie else: generateImages = False - rmg = loadRMGJob(inputFile, chemkinFile, speciesDict, generateImages=generateImages, useJava=java) + rmg = loadRMGJob(inputFile, chemkinFile, speciesDict, + generateImages=generateImages, useJava=java, checkDuplicates=checkDuplicates) if savePath is None: savePath = os.path.join(rmg.outputDirectory, 'flux') diff --git a/rmgpy/tools/loader.py b/rmgpy/tools/loader.py index d0c4e8a6d7..b72cc90b74 100644 --- a/rmgpy/tools/loader.py +++ b/rmgpy/tools/loader.py @@ -39,19 +39,23 @@ from rmgpy.solver.liquid import LiquidReactor from rmgpy.solver.base import TerminationConversion -def loadRMGJob(inputFile, chemkinFile=None, speciesDict=None, generateImages=True, useJava=False, useChemkinNames=False): +def loadRMGJob(inputFile, chemkinFile=None, speciesDict=None, generateImages=True, useJava=False, + useChemkinNames=False, checkDuplicates=True): if useJava: # The argument is an RMG-Java input file - rmg = loadRMGJavaJob(inputFile, chemkinFile, speciesDict, generateImages, useChemkinNames=useChemkinNames) + rmg = loadRMGJavaJob(inputFile, chemkinFile, speciesDict, generateImages, + useChemkinNames=useChemkinNames, checkDuplicates=checkDuplicates) else: # The argument is an RMG-Py input file - rmg = loadRMGPyJob(inputFile, chemkinFile, speciesDict, generateImages, useChemkinNames=useChemkinNames) + rmg = loadRMGPyJob(inputFile, chemkinFile, speciesDict, generateImages, + useChemkinNames=useChemkinNames, checkDuplicates=checkDuplicates) return rmg -def loadRMGPyJob(inputFile, chemkinFile=None, speciesDict=None, generateImages=True, useChemkinNames=False): +def loadRMGPyJob(inputFile, chemkinFile=None, speciesDict=None, generateImages=True, + useChemkinNames=False, checkDuplicates=True): """ Load the results of an RMG-Py job generated from the given `inputFile`. """ @@ -67,7 +71,8 @@ def loadRMGPyJob(inputFile, chemkinFile=None, speciesDict=None, generateImages=T chemkinFile = os.path.join(os.path.dirname(inputFile), 'chemkin', 'chem.inp') if not speciesDict: speciesDict = os.path.join(os.path.dirname(inputFile), 'chemkin', 'species_dictionary.txt') - speciesList, reactionList = loadChemkinFile(chemkinFile, speciesDict, useChemkinNames=useChemkinNames) + speciesList, reactionList = loadChemkinFile(chemkinFile, speciesDict, + useChemkinNames=useChemkinNames, checkDuplicates=checkDuplicates) # Map species in input file to corresponding species in Chemkin file speciesDict = {} @@ -121,7 +126,9 @@ def loadRMGPyJob(inputFile, chemkinFile=None, speciesDict=None, generateImages=T return rmg -def loadRMGJavaJob(inputFile, chemkinFile=None, speciesDict=None, useChemkinNames=False): + +def loadRMGJavaJob(inputFile, chemkinFile=None, speciesDict=None, generateImages=True, + useChemkinNames=False, checkDuplicates=True): """ Load the results of an RMG-Java job generated from the given `inputFile`. """ @@ -139,7 +146,8 @@ def loadRMGJavaJob(inputFile, chemkinFile=None, speciesDict=None, useChemkinName chemkinFile = os.path.join(os.path.dirname(inputFile), 'chemkin', 'chem.inp') if not speciesDict: speciesDict = os.path.join(os.path.dirname(inputFile), 'RMG_Dictionary.txt') - speciesList, reactionList = loadChemkinFile(chemkinFile, speciesDict, useChemkinNames=useChemkinNames) + speciesList, reactionList = loadChemkinFile(chemkinFile, speciesDict, + useChemkinNames=useChemkinNames, checkDuplicates=checkDuplicates) # Bath gas species don't appear in RMG-Java species dictionary, so handle # those as a special case @@ -178,15 +186,16 @@ def loadRMGJavaJob(inputFile, chemkinFile=None, speciesDict=None, useChemkinName rmg.reactionModel.core.reactions = reactionList # RMG-Java doesn't generate species images, so draw them ourselves now - speciesPath = os.path.join(os.path.dirname(inputFile), 'species') - try: - os.mkdir(speciesPath) - except OSError: - pass - for species in speciesList: - - path = os.path.join(speciesPath + '/{0!s}.png'.format(species)) - species.molecule[0].draw(str(path)) + if generateImages: + speciesPath = os.path.join(os.path.dirname(inputFile), 'species') + try: + os.mkdir(speciesPath) + except OSError: + pass + for species in speciesList: + path = os.path.join(speciesPath + '/{0!s}.png'.format(species)) + if not os.path.exists(path): + species.molecule[0].draw(str(path)) return rmg diff --git a/rmgpy/tools/sensitivity.py b/rmgpy/tools/simulate.py similarity index 90% rename from rmgpy/tools/sensitivity.py rename to rmgpy/tools/simulate.py index ce5a6c7869..3da1b14d0a 100644 --- a/rmgpy/tools/sensitivity.py +++ b/rmgpy/tools/simulate.py @@ -40,7 +40,7 @@ from rmgpy.solver.liquid import LiquidReactor from rmgpy.kinetics.diffusionLimited import diffusionLimiter -def plotSensitivity(outputDirectory, reactionSystemIndex, sensitiveSpeciesList, number=10, fileformat='.png'): +def plot_sensitivity(outputDirectory, reactionSystemIndex, sensitiveSpeciesList, number=10, fileformat='.png'): """ A function for plotting the top reaction thermo sensitivities (the number is inputted as the variable `number`) in bar plot format. @@ -93,12 +93,11 @@ def simulate(rmg, diffusionLimited=True): else: logging.info('Conducting simulation of reaction system %s...' % (index+1)) - - if rmg.saveSimulationProfiles: - reactionSystem.attach(SimulationProfileWriter( - rmg.outputDirectory, index, rmg.reactionModel.core.species)) - reactionSystem.attach(SimulationProfilePlotter( - rmg.outputDirectory, index, rmg.reactionModel.core.species)) + + reactionSystem.attach(SimulationProfileWriter( + rmg.outputDirectory, index, rmg.reactionModel.core.species)) + reactionSystem.attach(SimulationProfilePlotter( + rmg.outputDirectory, index, rmg.reactionModel.core.species)) sensWorksheet = [] for spec in reactionSystem.sensitiveSpecies: @@ -137,19 +136,19 @@ def simulate(rmg, diffusionLimited=True): ) if reactionSystem.sensitiveSpecies: - plotSensitivity(rmg.outputDirectory, index, reactionSystem.sensitiveSpecies) + plot_sensitivity(rmg.outputDirectory, index, reactionSystem.sensitiveSpecies) -def runSensitivity(inputFile, chemkinFile, dictFile, diffusionLimited=True): +def run_simulation(inputFile, chemkinFile, dictFile, diffusionLimited=True, checkDuplicates=True): """ Runs a standalone simulation of RMG. Runs sensitivity analysis if sensitive species are given. diffusionLimited=True implies that if it is a liquid reactor diffusion limitations will be enforced otherwise they will not be in a liquid reactor """ - rmg = loadRMGJob(inputFile, chemkinFile, dictFile, generateImages=False) + rmg = loadRMGJob(inputFile, chemkinFile, dictFile, generateImages=False, checkDuplicates=checkDuplicates) start_time = time() - # conduct sensitivity simulation + # conduct simulation simulate(rmg,diffusionLimited) end_time = time() time_taken = end_time - start_time diff --git a/rmgpy/tools/testSensitivity.py b/rmgpy/tools/simulateTest.py similarity index 92% rename from rmgpy/tools/testSensitivity.py rename to rmgpy/tools/simulateTest.py index 7db8dbbf38..d43c535931 100644 --- a/rmgpy/tools/testSensitivity.py +++ b/rmgpy/tools/simulateTest.py @@ -33,19 +33,19 @@ import os.path import shutil -from rmgpy.tools.sensitivity import runSensitivity +from rmgpy.tools.simulate import run_simulation import rmgpy -class SensitivityTest(unittest.TestCase): +class SimulateTest(unittest.TestCase): def test_minimal(self): - folder = os.path.join(os.path.dirname(rmgpy.__file__), 'tools/data/sens/simple') + folder = os.path.join(os.path.dirname(rmgpy.__file__), 'tools/data/sim/simple') inputFile = os.path.join(folder, 'input.py') chemkinFile = os.path.join(folder, 'chem.inp') dictFile = os.path.join(folder, 'species_dictionary.txt') - runSensitivity(inputFile, chemkinFile, dictFile) + run_simulation(inputFile, chemkinFile, dictFile) simfile = os.path.join(folder, 'solver', 'simulation_1_13.csv') sensfile = os.path.join(folder, 'solver', 'sensitivity_1_SPC_1.csv') @@ -56,13 +56,13 @@ def test_minimal(self): shutil.rmtree(os.path.join(folder, 'solver')) def test_liquid(self): - folder = os.path.join(os.path.dirname(rmgpy.__file__), 'tools/data/sens/liquid') + folder = os.path.join(os.path.dirname(rmgpy.__file__), 'tools/data/sim/liquid') inputFile = os.path.join(folder, 'input.py') chemkinFile = os.path.join(folder, 'chem.inp') dictFile = os.path.join(folder, 'species_dictionary.txt') - runSensitivity(inputFile, chemkinFile, dictFile, diffusionLimited=False) + run_simulation(inputFile, chemkinFile, dictFile, diffusionLimited=False) simfile = os.path.join(folder, 'solver', 'simulation_1_28.csv') sensfile = os.path.join(folder, 'solver', 'sensitivity_1_SPC_1.csv') diff --git a/rmgpy/tools/uncertainty.py b/rmgpy/tools/uncertainty.py index 11c0887983..a7d583057a 100644 --- a/rmgpy/tools/uncertainty.py +++ b/rmgpy/tools/uncertainty.py @@ -586,7 +586,7 @@ def sensitivityAnalysis(self, initialMoleFractions, sensitiveSpecies, T, P, term from rmgpy.solver import SimpleReactor, TerminationTime from rmgpy.quantity import Quantity - from rmgpy.tools.sensitivity import plotSensitivity + from rmgpy.tools.simulate import plot_sensitivity from rmgpy.rmg.listener import SimulationProfileWriter, SimulationProfilePlotter from rmgpy.rmg.settings import ModelSettings, SimulatorSettings T = Quantity(T) @@ -629,7 +629,7 @@ def sensitivityAnalysis(self, initialMoleFractions, sensitiveSpecies, T, P, term ) - plotSensitivity(self.outputDirectory, reactionSystemIndex, reactionSystem.sensitiveSpecies, number=number, fileformat=fileformat) + plot_sensitivity(self.outputDirectory, reactionSystemIndex, reactionSystem.sensitiveSpecies, number=number, fileformat=fileformat) def localAnalysis(self, sensitiveSpecies, correlated=False, number=10, fileformat='.png'): """ diff --git a/rmgpy/version.py b/rmgpy/version.py index bf6459b8cd..69e7c1ef82 100644 --- a/rmgpy/version.py +++ b/rmgpy/version.py @@ -34,4 +34,4 @@ This value can be accessed via `rmgpy.__version__`. """ -__version__ = '2.1.7' \ No newline at end of file +__version__ = '2.1.8' \ No newline at end of file diff --git a/scripts/generateChemkinHTML.py b/scripts/generateChemkinHTML.py index 21dcc6b5b0..8b6dfc033f 100644 --- a/scripts/generateChemkinHTML.py +++ b/scripts/generateChemkinHTML.py @@ -52,7 +52,8 @@ def main(chemkin, dictionary, output, foreign): model = CoreEdgeReactionModel() - model.core.species, model.core.reactions = loadChemkinFile(chemkin, dictionary, readComments=not foreign) + model.core.species, model.core.reactions = loadChemkinFile(chemkin, dictionary, readComments=not foreign, + checkDuplicates=foreign) outputPath = os.path.join(output, 'output.html') speciesPath = os.path.join(output, 'species') if not os.path.isdir(speciesPath): diff --git a/scripts/generateFluxDiagram.py b/scripts/generateFluxDiagram.py index 4ff457128b..719e8d5eb2 100644 --- a/scripts/generateFluxDiagram.py +++ b/scripts/generateFluxDiagram.py @@ -50,11 +50,13 @@ def parse_arguments(): parser.add_argument('input', metavar='INPUT', type=str, help='RMG input file') parser.add_argument('chemkin', metavar='CHEMKIN', type=str, help='Chemkin file') parser.add_argument('dictionary', metavar='DICTIONARY', type=str, help='RMG dictionary file') - parser.add_argument('species', metavar='SPECIES', type=str, nargs='?', default=None, help='Path to species images') parser.add_argument('chemkinOutput', metavar='CHEMKIN_OUTPUT', type=str, nargs='?', default=None, help='Chemkin output file') parser.add_argument('--java', action='store_true', help='process RMG-Java model') parser.add_argument('--no-dlim', dest='dlim', action='store_false', help='Turn off diffusion-limited rates') + parser.add_argument('-s', '--species', metavar='DIR', type=str, help='Path to folder containing species images') + parser.add_argument('-f', '--foreign', dest='checkDuplicates', action='store_true', + help='Not an RMG generated Chemkin file (will be checked for duplicates)') parser.add_argument('-n', '--maxnode', metavar='N', type=int, help='Maximum number of nodes to show in diagram') parser.add_argument('-e', '--maxedge', metavar='N', type=int, help='Maximum number of edges to show in diagram') parser.add_argument('-c', '--conctol', metavar='TOL', type=float, help='Lowest fractional concentration to show') @@ -71,18 +73,19 @@ def parse_arguments(): chemkinOutput = os.path.abspath(args.chemkinOutput) if args.chemkinOutput is not None else '' useJava = args.java dflag = args.dlim + checkDuplicates = args.checkDuplicates keys = ('maximumNodeCount', 'maximumEdgeCount', 'concentrationTolerance', 'speciesRateTolerance', 'timeStep') vals = (args.maxnode, args.maxedge, args.conctol, args.ratetol, args.tstep) settings = {k: v for k, v in zip(keys, vals) if v is not None} - return inputFile, chemkinFile, dictFile, speciesPath, chemkinOutput, useJava, dflag, settings + return inputFile, chemkinFile, dictFile, speciesPath, chemkinOutput, useJava, dflag, checkDuplicates, settings def main(): - inputFile, chemkinFile, dictFile, speciesPath, chemkinOutput, useJava, dflag, settings = parse_arguments() + inputFile, chemkinFile, dictFile, speciesPath, chemkinOutput, useJava, dflag, checkDuplicates, settings = parse_arguments() createFluxDiagram(inputFile, chemkinFile, dictFile, speciesPath=speciesPath, java=useJava, settings=settings, - chemkinOutput=chemkinOutput, diffusionLimited=dflag) + chemkinOutput=chemkinOutput, diffusionLimited=dflag, checkDuplicates=checkDuplicates) if __name__ == '__main__': main() \ No newline at end of file diff --git a/scripts/simulate.py b/scripts/simulate.py index f8705ae22c..f497e82a5b 100644 --- a/scripts/simulate.py +++ b/scripts/simulate.py @@ -29,14 +29,14 @@ ############################################################################### """ -This script runs stand-alone simulation on an RMG job. This is effectively the -same script as sensitivity.py +This script runs a stand-alone simulation (including sensitivity analysis if +specified in the input file) on an RMG job. """ import os.path import argparse -from rmgpy.tools.sensitivity import runSensitivity +from rmgpy.tools.simulate import run_simulation ################################################################################ @@ -49,27 +49,26 @@ def parse_arguments(): help='Chemkin file') parser.add_argument('dictionary', metavar='DICTIONARY', type=str, nargs=1, help='RMG dictionary file') + parser.add_argument('--no-dlim', dest='dlim', action='store_false', + help='Turn off diffusion-limited rates for LiquidReactor') + parser.add_argument('-f', '--foreign', dest='checkDuplicates', action='store_true', + help='Not an RMG generated Chemkin file (will be checked for duplicates)') args = parser.parse_args() inputFile = os.path.abspath(args.input[0]) chemkinFile = os.path.abspath(args.chemkin[0]) dictFile = os.path.abspath(args.dictionary[0]) + dflag = args.dlim + checkDuplicates = args.checkDuplicates - return inputFile, chemkinFile, dictFile + return inputFile, chemkinFile, dictFile, dflag, checkDuplicates def main(): - # This might not work anymore because functions were modified for use with webserver + inputFile, chemkinFile, dictFile, dflag, checkDuplicates = parse_arguments() - inputFile, chemkinFile, dictFile = parse_arguments() - - runSensitivity(inputFile, chemkinFile, dictFile) + run_simulation(inputFile, chemkinFile, dictFile, diffusionLimited=dflag, checkDuplicates=checkDuplicates) ################################################################################ - if __name__ == '__main__': main() - - - - \ No newline at end of file diff --git a/scripts/thermoEstimator.py b/scripts/thermoEstimator.py index 2cef325a12..19e3ffabc9 100644 --- a/scripts/thermoEstimator.py +++ b/scripts/thermoEstimator.py @@ -42,10 +42,11 @@ from rmgpy.chemkin import saveChemkinFile, saveSpeciesDictionary from rmgpy.rmg.model import Species from rmgpy.thermo.thermoengine import submit +from rmgpy.data.thermo import ThermoLibrary ################################################################################ -def runThermoEstimator(inputFile): +def runThermoEstimator(inputFile, library_flag): """ Estimate thermo for a list of species using RMG and the settings chosen inside a thermo input file. """ @@ -67,16 +68,17 @@ def runThermoEstimator(inputFile): for species in rmg.initialSpecies: submit(species) - # library = ThermoLibrary(name='Thermo Estimation Library') - # for spc in rmg.initialSpecies: - # library.loadEntry( - # index = len(library.entries) + 1, - # label = species.label, - # molecule = species.molecule[0].toAdjacencyList(), - # thermo = species.getThermoData().toThermoData(), - # shortDesc = species.getThermoData().comment, - # ) - # library.save(os.path.join(rmg.outputDirectory,'ThermoLibrary.py')) + if library_flag: + library = ThermoLibrary(name='Thermo Estimation Library') + for species in rmg.initialSpecies: + library.loadEntry( + index = len(library.entries) + 1, + label = species.label, + molecule = species.molecule[0].toAdjacencyList(), + thermo = species.getThermoData().toThermoData(), + shortDesc = species.getThermoData().comment, + ) + library.save(os.path.join(rmg.outputDirectory,'ThermoLibrary.py')) # Save the thermo data to chemkin format output files and dictionary, with no reactions @@ -92,8 +94,10 @@ def runThermoEstimator(inputFile): parser = argparse.ArgumentParser() parser.add_argument('input', metavar='INPUT', type=str, nargs=1, help='Thermo input file') + parser.add_argument('-l', '--library', action='store_true', help='generate RMG thermo library') + args = parser.parse_args() inputFile = os.path.abspath(args.input[0]) - runThermoEstimator(inputFile) \ No newline at end of file + runThermoEstimator(inputFile, args.library) diff --git a/setup.py b/setup.py index 75df912420..58cf01b0d2 100644 --- a/setup.py +++ b/setup.py @@ -109,6 +109,7 @@ def getMainExtensionModules(): Extension('rmgpy.quantity', ['rmgpy/quantity.py'], include_dirs=['.']), Extension('rmgpy.reaction', ['rmgpy/reaction.py'], include_dirs=['.']), Extension('rmgpy.species', ['rmgpy/species.py'], include_dirs=['.']), + Extension('rmgpy.chemkin', ['rmgpy/chemkin.pyx'], include_dirs=['.']), ] def getSolverExtensionModules(): @@ -185,8 +186,19 @@ def getCanthermExtensionModules(): if os.path.splitext(source)[1] == '.pyx': ext_modules.append(module) -scripts=['cantherm.py', 'rmg.py', 'scripts/diffModels.py', 'scripts/generateFluxDiagram.py', - 'scripts/generateReactions.py', 'scripts/mergeModels.py','scripts/sensitivity.py', 'scripts/thermoEstimator.py', +scripts=['cantherm.py', + 'rmg.py', + 'scripts/checkModels.py', + 'scripts/convertFAME.py', + 'scripts/diffModels.py', + 'scripts/generateChemkinHTML.py', + 'scripts/generateFluxDiagram.py', + 'scripts/generateReactions.py', + 'scripts/machineWriteDatabase.py', + 'scripts/mergeModels.py', + 'scripts/simulate.py', + 'scripts/standardizeModelSpeciesNames.py', + 'scripts/thermoEstimator.py', 'testing/databaseTest.py'] modules = [] diff --git a/update_headers.py b/update_headers.py deleted file mode 100644 index ad01215e15..0000000000 --- a/update_headers.py +++ /dev/null @@ -1,112 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- - -""" -This script automatically updates the header for *.py, *.pyx, and *.pxd -files in RMG-Py/rmgpy, RMG-Py/scripts, and the root RMG-Py directory. - -Other directories (e.g., examples, documentation, etc.) are ignored -along with directories containing test data. - -The headers are automatically generated based on the LICENSE.txt file. -Typical usage would involve running this script after updating the copyright -date in the license file. - -Because this script makes assumptions regarding the contents at the -start of each file, be sure to double-check the results to make sure -important lines aren't accidentally overwritten. -""" - -import os - -shebang = """#!/usr/bin/env python -# -*- coding: utf-8 -*- - -""" - -header = """############################################################################### -# # -# RMG - Reaction Mechanism Generator # -# # -""" - -with open('LICENSE.txt', 'r') as f: - for line in f: - line = line.strip() - newline = '# {0:<75} #\n'.format(line) - header += newline - -header += """# # -############################################################################### -""" - -print header - -def replace_header(oldfile): - newfile = os.path.join('tmp', 'tempfile') - - ext = os.path.splitext(oldfile)[1] - if ext == '.py': - py_file = True - elif ext == '.pyx' or ext == '.pxd': - py_file = False - else: - raise Exception('Unexpected file type: {0}'.format(oldfile)) - - with open(oldfile, 'r') as old, open(newfile, 'w+') as new: - - # Write shebang and encoding for python files only - if py_file: - new.write(shebang) - - # Read old file and copy over contents - found_bar = False - first_line = True - start = False - for i, line in enumerate(old): - if i == 0 and line[0] != '#': - # Assume there's no header, so start copying lines right away - new.write(header) - start = True - if start: - if first_line and line.strip() != '': - new.write('\n') - first_line = False - new.write(line) - elif line.startswith('# cython:'): - # Copy over any cython directives - new.write(line) - new.write('\n') - elif line.startswith('##########'): - if found_bar: - # We've reached the end of the license, so start copying lines starting with the next line - new.write(header) - start = True - else: - found_bar = True - - # Replace old file with new file - os.rename(newfile, oldfile) - -# Create temporary directory for writing files -if not os.path.exists('tmp'): - os.makedirs('tmp') - -# Compile list of files to modify -filelist = ['rmg.py', 'cantherm.py', 'setup.py'] - -root_dirs = ['rmgpy', 'scripts'] -for root_dir in root_dirs: - for root, dirs, files in os.walk(root_dir): - if 'test_data' in root or 'files' in root or '/tools/data' in root or '/cantherm/data' in root: - continue - print root - for f in files: - ext = os.path.splitext(f)[1] - if ext in ['.py', '.pyx', '.pxd']: - filelist.append(os.path.join(root, f)) - -for f in filelist: - print 'Updating {0} ...'.format(f) - replace_header(f) - diff --git a/utilities.py b/utilities.py new file mode 100644 index 0000000000..cb5f14b7b0 --- /dev/null +++ b/utilities.py @@ -0,0 +1,292 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +############################################################################### +# # +# RMG - Reaction Mechanism Generator # +# # +# Copyright (c) 2002-2018 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 # +# copy of this software and associated documentation files (the 'Software'), # +# to deal in the Software without restriction, including without limitation # +# the rights to use, copy, modify, merge, publish, distribute, sublicense, # +# and/or sell copies of the Software, and to permit persons to whom the # +# Software is furnished to do so, subject to the following conditions: # +# # +# The above copyright notice and this permission notice shall be included in # +# all copies or substantial portions of the Software. # +# # +# THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR # +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE # +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER # +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING # +# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER # +# DEALINGS IN THE SOFTWARE. # +# # +############################################################################### + +""" +This module contains utilities related to compilation and code maintenance. +""" + +import argparse +import os +import platform +import re +import subprocess + + +def check_dependencies(): + """ + Checks for and locates major dependencies that RMG requires. + """ + missing = False + + print '\nChecking vital dependencies...\n' + print '{0:<15}{1:<15}{2}'.format('Package', 'Version', 'Location') + + # Check for symmetry + try: + result = subprocess.check_output('symmetry -h', stderr=subprocess.STDOUT, shell=True) + except subprocess.CalledProcessError: + print '{0:<30}{1}'.format('symmetry', 'Not found. Please install in order to use QM.') + missing = True + else: + match = re.search(r'\$Revision: (\S*) \$', result) + version = match.group(1) + if platform.system() == 'Windows': + location = subprocess.check_output('where symmetry', shell=True) + else: + location = subprocess.check_output('which symmetry', shell=True) + + print '{0:<15}{1:<15}{2}'.format('symmetry', version, location.strip()) + + # Check for RDKit + try: + import rdkit + from rdkit import Chem + except ImportError: + print '{0:<30}{1}'.format('RDKit', + 'Not found. Please install RDKit version 2015.03.1 or later with InChI support.') + missing = True + else: + try: + version = rdkit.__version__ + except AttributeError: + version = False + location = rdkit.__file__ + inchi = Chem.inchi.INCHI_AVAILABLE + + if version: + print '{0:<15}{1:<15}{2}'.format('RDKit', version, location) + if not inchi: + print 'RDKit installed without InChI Support. Please install with InChI.' + missing = True + else: + print 'RDKit version out of date, please install RDKit version 2015.03.1 or later with InChI support.' + missing = True + + # Check for OpenBabel + try: + import openbabel + except ImportError: + print '{0:<30}{1}'.format('OpenBabel', + 'Not found. Necessary for SMILES/InChI functionality for nitrogen compounds.') + missing = True + else: + location = openbabel.__file__ + print '{0:<30}{1}'.format('OpenBabel', location) + + # Check for lpsolve + 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) + + # Check for pyrdl + try: + import py_rdl + except ImportError: + print '{0:<30}{1}'.format('pyrdl', 'Not found. Necessary for ring perception algorithms.') + missing = True + else: + location = py_rdl.__file__ + print '{0:<30}{1}'.format('pyrdl', location) + + if missing: + print """ +There are missing dependencies as listed above. Please install them before proceeding. + +Using Anaconda, these dependencies can be individually installed from the RMG channel as follows: + + conda install -c rmg [package name] + +You can alternatively update your environment and install all missing dependencies as follows: + + conda env update -f environment_[linux/mac/windows].yml # Choose the correct file for your OS + +Be sure to activate your conda environment (rmg_env by default) before installing or updating. +""" + else: + print """ +Everything was found :) +""" + + +def clean(subdirectory=''): + """ + Removes files generated during compilation. + + For *nix systems, remove `.so` files and `.pyc` files. + For Windows, remove `.pyd` files and `.pyc` files. + """ + if platform.system() == 'Windows': + extensions = ['.pyd', '.pyc'] + else: + extensions = ['.so', '.pyc'] + + directory = os.path.join(os.path.dirname(os.path.abspath(__file__)), subdirectory) + + for root, dirs, files in os.walk(directory): + for f in files: + ext = os.path.splitext(f)[1] + if ext in extensions: + os.remove(os.path.join(root, f)) + + +def update_headers(): + """ + Automatically update the headers for *.py, *.pyx, and *.pxd files + in RMG-Py/rmgpy, RMG-Py/scripts, and the root RMG-Py directory. + + Other directories (e.g., examples, documentation, etc.) are ignored + along with directories containing test data. + + The headers are automatically generated based on the LICENSE.txt file. + Typical usage would involve running this script after updating the + copyright date in the license file. + + Because this script makes assumptions regarding the contents at the + 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 +# -*- coding: utf-8 -*- + +""" + + header = """############################################################################### +# # +# RMG - Reaction Mechanism Generator # +# # +""" + + with open('LICENSE.txt', 'r') as f: + for line in f: + line = line.strip() + newline = '# {0:<75} #\n'.format(line) + header += newline + + header += """# # +############################################################################### +""" + + print header + + def replace_header(oldfile): + newfile = os.path.join('tmp', 'tempfile') + + ext = os.path.splitext(oldfile)[1] + if ext == '.py': + py_file = True + elif ext == '.pyx' or ext == '.pxd': + py_file = False + else: + raise Exception('Unexpected file type: {0}'.format(oldfile)) + + with open(oldfile, 'r') as old, open(newfile, 'w+') as new: + + # Write shebang and encoding for python files only + if py_file: + new.write(shebang) + + # Read old file and copy over contents + found_bar = False + first_line = True + start = False + for i, line in enumerate(old): + if i == 0 and line[0] != '#': + # Assume there's no header, so start copying lines right away + new.write(header) + start = True + if start: + if first_line and line.strip() != '': + new.write('\n') + first_line = False + new.write(line) + elif line.startswith('# cython:'): + # Copy over any cython directives + new.write(line) + new.write('\n') + elif line.startswith('##########'): + if found_bar: + # We've reached the end of the license, so start copying lines starting with the next line + new.write(header) + start = True + else: + found_bar = True + + # Replace old file with new file + os.rename(newfile, oldfile) + + # Create temporary directory for writing files + if not os.path.exists('tmp'): + os.makedirs('tmp') + + # Compile list of files to modify + filelist = ['rmg.py', 'cantherm.py', 'setup.py'] + + root_dirs = ['rmgpy', 'scripts'] + for root_dir in root_dirs: + for root, dirs, files in os.walk(root_dir): + if 'test_data' in root or 'files' in root or '/tools/data' in root or '/cantherm/data' in root: + continue + print root + for f in files: + ext = os.path.splitext(f)[1] + if ext in ['.py', '.pyx', '.pxd']: + filelist.append(os.path.join(root, f)) + + for f in filelist: + print 'Updating {0} ...'.format(f) + replace_header(f) + + # Remove temporary directory + os.rmdir('tmp') + + +if __name__ == '__main__': + parser = argparse.ArgumentParser(description='RMG Code Utilities') + + parser.add_argument('command', metavar='COMMAND', type=str, + choices=['check-dependencies', 'clean', 'clean-solver', 'update-headers'], + help='command to execute') + + args = parser.parse_args() + + if args.command == 'check-dependencies': + check_dependencies() + elif args.command == 'clean': + clean() + elif args.command == 'clean-solver': + clean('rmgpy/solver') + elif args.command == 'update-headers': + update_headers()