Surface auto tree gen support with metal attrs

.github/workflows/CI.yml

@@ -12,7 +12,7 @@ jobs:
     strategy:
       max-parallel: 5
     env: # update this if needed to match a pull request on the RMG-database
-      RMG_DATABASE_BRANCH: main
+      RMG_DATABASE_BRANCH: cathub
     defaults:
       run:
         shell: bash -l {0}

ipython/estimate_surface_A_factors.ipynb

@@ -0,0 +1,154 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "4a175662",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from rmgpy.data.rmg import RMGDatabase\n",
+    "from rmgpy.kinetics import MultiArrhenius\n",
+    "from rmgpy import settings\n",
+    "\n",
+    "from surface_A_factor import *\n",
+    "from IPython.display import display\n",
+    "import numpy as np\n",
+    "import os"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "id": "a08414a4",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "['Surface/cathub/Pd',\n",
+       " 'Surface/cathub/Ni',\n",
+       " 'Surface/cathub/Ir',\n",
+       " 'Surface/cathub/Cu',\n",
+       " 'Surface/cathub/Co',\n",
+       " 'Surface/cathub/Fe',\n",
+       " 'Surface/cathub/Ru',\n",
+       " 'Surface/cathub/Rh',\n",
+       " 'Surface/cathub/Pt',\n",
+       " 'Surface/cathub/Ag',\n",
+       " 'Surface/cathub/Au',\n",
+       " 'Surface/cathub/W']"
+      ]
+     },
+     "execution_count": 6,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "libs = []\n",
+    "p = 'Surface/cathub/'\n",
+    "for e in os.listdir(os.path.join(settings['database.directory'],'kinetics','libraries',p)):\n",
+    "    libs.append(os.path.join(p,e))\n",
+    "libs"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "id": "84ef9667",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "WARNING:root:An instance of RMGDatabase already exists. Re-initializing it.\n"
+     ]
+    }
+   ],
+   "source": [
+    "database = RMGDatabase()\n",
+    "database.load(\n",
+    "    settings['database.directory'],\n",
+    "    thermo_libraries=[\n",
+    "     'primaryThermoLibrary','DFT_QCI_thermo','surfaceThermoPt111','thermo_DFT_CCSDTF12_BAC',\n",
+    "        'CHON_G4', 'NISTThermoLibrary','BurcatNS','JetSurF2.0',\n",
+    "    ],# can add more\n",
+    "    transport_libraries=[],\n",
+    "    reaction_libraries=libs,\n",
+    "    seed_mechanisms=[],\n",
+    "    kinetics_families=\"None\",\n",
+    "    kinetics_depositories=[],\n",
+    "    statmech_libraries=[],\n",
+    "    depository=[],\n",
+    "    solvation=False,\n",
+    "    surface=True,\n",
+    "    testing=False,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "id": "fa5d6520",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "processing Surface/cathub/Pd...\n",
+      "processing Surface/cathub/Ni...\n",
+      "processing Surface/cathub/Ir...\n",
+      "processing Surface/cathub/Cu...\n",
+      "processing Surface/cathub/Co...\n",
+      "processing Surface/cathub/Fe...\n",
+      "processing Surface/cathub/Ru...\n",
+      "processing Surface/cathub/Rh...\n",
+      "processing Surface/cathub/Pt...\n",
+      "processing Surface/cathub/Ag...\n",
+      "processing Surface/cathub/Au...\n",
+      "processing Surface/cathub/W...\n"
+     ]
+    }
+   ],
+   "source": [
+    "for name,library in database.kinetics.libraries.items():\n",
+    "    print(f'processing {name}...')\n",
+    "    for i,entry in library.entries.items():\n",
+    "        A,comment = estimate_surface_A(entry.item)\n",
+    "        entry.long_desc += '\\n'\n",
+    "        entry.long_desc += comment\n",
+    "        if isinstance(entry.data,MultiArrhenius):\n",
+    "            for arr in entry.data.arrhenius:\n",
+    "                arr.A = A\n",
+    "        else:\n",
+    "            entry.data.A = A\n",
+    "    path = os.path.join(settings['database.directory'],'kinetics','libraries',name,'reactions.py')\n",
+    "    library.save(path)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.10"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

ipython/surface_A_factor.py

@@ -0,0 +1,172 @@
+#!/usr/bin/env python3
+
+from rmgpy.molecule import Molecule
+from rmgpy.species import Species
+from rmgpy.constants import kB,h,R
+from rmgpy.thermo.thermoengine import submit
+from rmgpy.kinetics import get_rate_coefficient_units_from_reaction_order
+
+import logging
+import numpy as np
+
+def get_A_units(rxn):
+
+    try:
+        surf_reacts = [spcs for spcs in rxn.reactants if spcs.contains_surface_site()]
+    except IndexError:
+        surf_prods = []
+        logging.warning(f"Species do not have an rmgpy.molecule.Molecule "  
+                        "Cannot determine phases of species. We will assume gas"
+                        )
+    n_surf = len(surf_reacts)
+    n_gas = len(rxn.reactants) - len(surf_reacts)
+    return get_rate_coefficient_units_from_reaction_order(n_gas, n_surf)
+
+def get_surface_reaction_stoich(rxn):
+
+    n_gas = 0
+    n_ads = 0
+    n_x = 0
+
+    for s in rxn.reactants:
+        if s.is_surface_site():
+            n_x -= 1
+        elif s.contains_surface_site():
+            n_ads -= 1
+        else:
+            n_gas -= 1
+    for s in rxn.products:
+        if s.is_surface_site():
+            n_x += 1
+        elif s.contains_surface_site():
+            n_ads += 1
+        else:
+            n_gas += 1
+
+    return (n_gas,n_ads,n_x)
+
+def get_surface_reaction_type(rxn):
+
+    n_gas, n_ads, n_x = get_surface_reaction_stoich(rxn)
+
+    if (n_gas, n_ads, n_x) == (1,-1,1):
+        return "desorption"
+    elif (n_gas, n_ads, n_x) == (-1,1,-1):
+        return "adsorption"
+    elif (n_gas, n_ads, n_x) == (0,1,-1):
+        return "dissociation"
+    elif (n_gas, n_ads, n_x) == (0,-1,1):
+        return "association"
+    else:
+        return "unknown"
+
+def estimate_surface_desorption_A(adsorbate: Species):
+
+    Ar_mw = 39.87750372310267 # amu
+    dummy_species = get_dummy_species(adsorbate)
+    S = dummy_species.get_entropy(298.)
+    mw = dummy_species.molecular_weight.value
+
+    A = kB*298./h
+    A *= np.exp(0.3*S/R+3.3-(18.6+np.log((mw/Ar_mw)**(5./2.)))/3)
+    commment = f"A factor estimated from gas-phase smiles {dummy_species.smiles} from "\
+    f"{dummy_species.thermo.comment} and S298={S/4.184:.2f} cal/mol/K"
+    return A,commment
+
+def get_dummy_species(adsorbate: Species):
+
+    dummy_molecules = adsorbate.molecule[0].get_desorbed_molecules()
+    for mol in dummy_molecules:
+        mol.clear_labeled_atoms()
+
+    gas_phase_species_from_libraries = []
+    gas_phase_species_estimates = []
+    for dummy_molecule in dummy_molecules:
+        dummy_species = Species()
+        dummy_species.molecule = [dummy_molecule]
+        dummy_species.generate_resonance_structures()
+        submit(dummy_species)
+        if dummy_species.thermo.label:
+            gas_phase_species_from_libraries.append(dummy_species)
+        else:
+            gas_phase_species_estimates.append(dummy_species)
+
+    # define the comparison function to find the lowest energy
+    def lowest_energy(species):
+        if hasattr(species.thermo, 'H298'):
+            return species.thermo.H298.value_si
+        else:
+            return species.thermo.get_enthalpy(298.0)
+
+    if gas_phase_species_from_libraries:
+        species = min(gas_phase_species_from_libraries, key=lowest_energy)
+    else:
+        species = min(gas_phase_species_estimates, key=lowest_energy)
+
+    return species
+
+def estimate_surface_dissociation_A(adsorbate: Species):
+
+    A, comment = estimate_surface_desorption_A(adsorbate)
+    return 0.001 * A, comment
+
+def get_adsorbate(rxn, on='reactants'):
+
+    if on == 'reactants':
+        for s in rxn.reactants:
+            if s.contains_surface_site() and not s.is_surface_site():
+                return s
+    elif on == 'products':
+        for s in rxn.products:
+            if s.contains_surface_site() and not s.is_surface_site():
+                return s
+
+def estimate_surface_A(rxn):
+
+    reaction_type = get_surface_reaction_type(rxn)
+    comment = "A factor estimation:\n"
+
+    if reaction_type == 'desorption':
+        comment += f"A factor estimate for {reaction_type}\n"
+        adsorbate = get_adsorbate(rxn,'reactants')
+        A, _comment = estimate_surface_desorption_A(adsorbate)
+        comment += _comment
+
+    elif reaction_type == 'dissociation':
+        comment += f"A factor estimate for {reaction_type}\n"
+        adsorbate = get_adsorbate(rxn,'reactants')
+        A, _comment = estimate_surface_dissociation_A(adsorbate)
+        comment += _comment
+
+    elif reaction_type == 'association':
+        comment += f"A factor estimate for {reaction_type}\n"
+        adsorbate = get_adsorbate(rxn,'products')
+        A, _comment = estimate_surface_dissociation_A(adsorbate)
+        comment += _comment
+        if not adsorbate.thermo:
+            submit(adsorbate)
+        deltaS = adsorbate.thermo.get_entropy(298.)
+        assert len(rxn.reactants) == 2
+        for s in rxn.reactants:
+            if s.contains_surface_site():
+                if not s.thermo:
+                    submit(s)
+                deltaS -= s.thermo.get_entropy(298.)
+        A *= np.exp(deltaS/R)
+
+    else:
+        A = kB*298./h
+        comment = "Could not determine reaction type "\
+        f"estimating A = kb/298/h = {A:.2e}"
+
+    units = get_A_units(rxn)
+    surface_site_density=2.483e-05 # Pt111 in metal DB
+    if units in ('m^2/(mol*s)','m^5/(mol^2*s)'):
+        A /= surface_site_density
+        comment += '\nA/=2.483e-5 mol/m^2 (Pt111 site density)'
+    elif units == 'm^4/(mol^2*s)':
+        A /= surface_site_density
+        A /= surface_site_density
+        comment += '\nA/=(2.483e-5 mol/m^2)^2 (Pt111 site density)'
+
+    return (A,units), comment