From f010e4060922576c56a98517205fda74f390c52f Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 22 Jul 2021 11:52:46 -0400 Subject: [PATCH 01/36] added `get_metal_label` entry method --- rmgpy/data/base.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/rmgpy/data/base.py b/rmgpy/data/base.py index 2b617a6b27..5850bd3bad 100644 --- a/rmgpy/data/base.py +++ b/rmgpy/data/base.py @@ -126,6 +126,16 @@ def __str__(self): def __repr__(self): return ''.format(self.index, self.label) + def get_metal_label(self): + """ + retrieve the metal label (str) for the entry + """ + if self.facet is None: + metal = self.metal # could be None + else: + metal = self.metal + self.facet + return metal + def get_all_descendants(self): """ retrieve all the descendants of entry From f9a43e3a8558d280f4352712acb3c2407d88f1b2 Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 22 Jul 2021 11:58:40 -0400 Subject: [PATCH 02/36] use `get_rate_coefficient_units_from_reaction_order` to get A units for ArrBM this is needed so we can get the correct units for surface reactions --- rmgpy/kinetics/arrhenius.pyx | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/rmgpy/kinetics/arrhenius.pyx b/rmgpy/kinetics/arrhenius.pyx index 244318a61d..98377f42da 100644 --- a/rmgpy/kinetics/arrhenius.pyx +++ b/rmgpy/kinetics/arrhenius.pyx @@ -36,6 +36,7 @@ cimport rmgpy.constants as constants import rmgpy.quantity as quantity from rmgpy.exceptions import KineticsError from rmgpy.kinetics.uncertainties import rank_accuracy_map +from rmgpy.kinetics import get_rate_coefficient_units_from_reaction_order from rmgpy.molecule.molecule import Bond # Prior to numpy 1.14, `numpy.linalg.lstsq` does not accept None as a value @@ -669,9 +670,11 @@ cdef class ArrheniusBM(KineticsModel): self.comment = 'Fitted to {0} reactions at temperatures: {1}'.format(len(rxns), Ts) # fill in parameters - A_units = ['', 's^-1', 'm^3/(mol*s)', 'm^6/(mol^2*s)'] - order = len(rxns[0].reactants) - self.A = (A, A_units[order]) + surf_reacts = [spcs for spcs in rxns[0].reactants if spcs.contains_surface_site()] + n_surf = len(surf_reacts) + n_gas = len(rxns[0].reactants) - len(surf_reacts) + A_units = get_rate_coefficient_units_from_reaction_order(n_gas, n_surf) + self.A = (A, A_units) self.n = n self.w0 = (w0, 'J/mol') From dae5b644bda98f1c86b01462fa2e3402efa5cb68 Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 22 Jul 2021 12:00:09 -0400 Subject: [PATCH 03/36] `get_w0` with metals and vdw bonds --- rmgpy/kinetics/arrhenius.pyx | 29 +++++++++++++++++++++-------- 1 file changed, 21 insertions(+), 8 deletions(-) diff --git a/rmgpy/kinetics/arrhenius.pyx b/rmgpy/kinetics/arrhenius.pyx index 98377f42da..f0229c1213 100644 --- a/rmgpy/kinetics/arrhenius.pyx +++ b/rmgpy/kinetics/arrhenius.pyx @@ -1109,6 +1109,8 @@ def get_w0(actions, rxn): recipe = actions + metal = rxn.get_metal_label() + wb = 0.0 wf = 0.0 for act in recipe: @@ -1123,14 +1125,24 @@ def get_w0(actions, rxn): atom2 = a_dict[act[3]] if act[0] == 'BREAK_BOND': - bd = mol.get_bond(atom1, atom2) - wb += bd.get_bde() + bd = mol.get_bond(a_dict[act[1]], a_dict[act[3]]) + wb += bd.get_bde(metal=metal) elif act[0] == 'FORM_BOND': - bd = Bond(atom1, atom2, act[2]) - wf += bd.get_bde() + bd = Bond(a_dict[act[1]], a_dict[act[3]], act[2]) + wf += bd.get_bde(metal=metal) elif act[0] == 'CHANGE_BOND': - bd1 = mol.get_bond(atom1, atom2) - + if not mol.has_bond(a_dict[act[1]], a_dict[act[3]]): # we dont have a bond + is_vdW_bond = False + for atom in (a_dict[act[1]], a_dict[act[3]]): + if atom.is_surface_site(): + is_vdW_bond = True + break + if not is_vdW_bond: # no surface site, so no vdW bond + raise ('Attempted to change a nonexistent bond.') + else: # we found a surface site, so we will make vdw bond + bd1 = Bond(a_dict[act[1]], a_dict[act[3]], order=0) + else: # we have a bond + bd1 = mol.get_bond(a_dict[act[1]], a_dict[act[3]]) if act[2] + bd1.order == 0.5: mol2 = None for r in rxn.products: @@ -1154,8 +1166,9 @@ def get_w0(actions, rxn): if bd2.order == 0: bd2_bde = 0.0 else: - bd2_bde = bd2.get_bde() - bde_diff = bd2_bde - bd1.get_bde() + bd2_bde = bd2.get_bde(metal=metal) + + bde_diff = bd2_bde - bd1.get_bde(metal=metal) if bde_diff > 0: wf += abs(bde_diff) else: From 606e3c025a7967f5997dcba48954270cdaa95875 Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 22 Jul 2021 12:03:37 -0400 Subject: [PATCH 04/36] added SurfaceArrheniusBM kinetics class --- rmgpy/data/kinetics/database.py | 8 ++-- rmgpy/data/kinetics/family.py | 2 +- rmgpy/kinetics/__init__.py | 2 +- rmgpy/kinetics/surface.pxd | 5 ++- rmgpy/kinetics/surface.pyx | 73 +++++++++++++++++++++++++++++++++ rmgpy/reaction.py | 2 +- 6 files changed, 85 insertions(+), 7 deletions(-) diff --git a/rmgpy/data/kinetics/database.py b/rmgpy/data/kinetics/database.py index d37e0de432..7525f5ed62 100644 --- a/rmgpy/data/kinetics/database.py +++ b/rmgpy/data/kinetics/database.py @@ -44,7 +44,8 @@ from rmgpy.kinetics import Arrhenius, ArrheniusEP, ThirdBody, Lindemann, Troe, \ PDepArrhenius, MultiArrhenius, MultiPDepArrhenius, \ Chebyshev, KineticsData, StickingCoefficient, \ - StickingCoefficientBEP, SurfaceArrhenius, SurfaceArrheniusBEP, ArrheniusBM + StickingCoefficientBEP, SurfaceArrhenius, SurfaceArrheniusBEP, \ + ArrheniusBM, SurfaceArrheniusBM from rmgpy.molecule import Molecule, Group from rmgpy.reaction import Reaction, same_species_lists from rmgpy.species import Species @@ -79,7 +80,8 @@ def __init__(self): 'SurfaceArrhenius': SurfaceArrhenius, 'SurfaceArrheniusBEP': SurfaceArrheniusBEP, 'R': constants.R, - 'ArrheniusBM': ArrheniusBM + 'ArrheniusBM': ArrheniusBM, + 'SurfaceArrheniusBM': SurfaceArrheniusBM } self.global_context = {} @@ -813,6 +815,6 @@ def reconstruct_kinetics_from_source(self, reaction, source, fix_barrier_height= else: h298 = rxn_copy.get_enthalpy_of_reaction(298) - if isinstance(kinetics, (ArrheniusEP, ArrheniusBM)): + if isinstance(kinetics, (ArrheniusEP, ArrheniusBM, SurfaceArrheniusBM)): kinetics = kinetics.to_arrhenius(h298) return kinetics diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 29794d7875..0a7b14ca2e 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -55,7 +55,7 @@ from rmgpy.exceptions import ActionError, DatabaseError, InvalidActionError, KekulizationError, KineticsError, \ ForbiddenStructureException, UndeterminableKineticsError from rmgpy.kinetics import Arrhenius, SurfaceArrhenius, SurfaceArrheniusBEP, StickingCoefficient, \ - StickingCoefficientBEP, ArrheniusBM + StickingCoefficientBEP, ArrheniusBM, SurfaceArrheniusBM from rmgpy.kinetics.uncertainties import RateUncertainty, rank_accuracy_map from rmgpy.molecule import Bond, GroupBond, Group, Molecule from rmgpy.molecule.atomtype import ATOMTYPES diff --git a/rmgpy/kinetics/__init__.py b/rmgpy/kinetics/__init__.py index a9dc32e49b..10bc31ef48 100644 --- a/rmgpy/kinetics/__init__.py +++ b/rmgpy/kinetics/__init__.py @@ -35,4 +35,4 @@ from rmgpy.kinetics.kineticsdata import KineticsData, PDepKineticsData from rmgpy.kinetics.tunneling import Wigner, Eckart from rmgpy.kinetics.surface import SurfaceArrhenius, SurfaceArrheniusBEP, \ - StickingCoefficient, StickingCoefficientBEP + StickingCoefficient, StickingCoefficientBEP, SurfaceArrheniusBM diff --git a/rmgpy/kinetics/surface.pxd b/rmgpy/kinetics/surface.pxd index 8e37ac54b7..bbd87ca5b9 100644 --- a/rmgpy/kinetics/surface.pxd +++ b/rmgpy/kinetics/surface.pxd @@ -28,7 +28,7 @@ cimport numpy as np from rmgpy.kinetics.model cimport KineticsModel -from rmgpy.kinetics.arrhenius cimport Arrhenius, ArrheniusEP +from rmgpy.kinetics.arrhenius cimport Arrhenius, ArrheniusEP, ArrheniusBM from rmgpy.quantity cimport ScalarQuantity, ArrayQuantity ################################################################################ @@ -72,4 +72,7 @@ cdef class SurfaceArrhenius(Arrhenius): cdef class SurfaceArrheniusBEP(ArrheniusEP): cdef public dict _coverage_dependence pass +################################################################################ +cdef class SurfaceArrheniusBM(ArrheniusBM): + pass diff --git a/rmgpy/kinetics/surface.pyx b/rmgpy/kinetics/surface.pyx index 2c108e9204..3bf48cd4f0 100644 --- a/rmgpy/kinetics/surface.pyx +++ b/rmgpy/kinetics/surface.pyx @@ -625,3 +625,76 @@ cdef class SurfaceArrheniusBEP(ArrheniusEP): coverage_dependence=self.coverage_dependence, comment=self.comment, ) + +################################################################################ + +cdef class SurfaceArrheniusBM(ArrheniusBM): + """ + A kinetics model based on the (modified) Arrhenius equation, using the + Blowers-Masel equation to determine the activation energy. + Based on Blowers and Masel's 2000 paper Engineering Approximations for Activation + Energies in Hydrogen Transfer Reactions. + The attributes are: + + =============== ============================================================= + Attribute Description + =============== ============================================================= + `A` The preexponential factor + `n` The temperature exponent + `w0` The average of the bond dissociation energies of the bond formed and the bond broken + `E0` The activation energy for a thermoneutral reaction + `Tmin` The minimum temperature at which the model is valid, or zero if unknown or undefined + `Tmax` The maximum temperature at which the model is valid, or zero if unknown or undefined + `Pmin` The minimum pressure at which the model is valid, or zero if unknown or undefined + `Pmax` The maximum pressure at which the model is valid, or zero if unknown or undefined + `comment` Information about the model (e.g. its source) + =============== ============================================================= + + """ + + property A: + """The preexponential factor. + + This is the only thing different from a ArrheniusBM class""" + def __get__(self): + return self._A + def __set__(self, value): + self._A = quantity.SurfaceRateCoefficient(value) + + def __repr__(self): + """ + Return a string representation that can be used to reconstruct the + ArrheniusBM object. + """ + string = 'SurfaceArrheniusBM(A={0!r}, n={1!r}, w0={2!r}, E0={3!r}'.format(self.A, self.n, self.w0, self.E0) + if self.Tmin is not None: string += ', Tmin={0!r}'.format(self.Tmin) + if self.Tmax is not None: string += ', Tmax={0!r}'.format(self.Tmax) + if self.Pmin is not None: string += ', Pmin={0!r}'.format(self.Pmin) + if self.Pmax is not None: string += ', Pmax={0!r}'.format(self.Pmax) + if self.uncertainty is not None: string += ', uncertainty={0!r}'.format(self.uncertainty) + if self.comment != '': string += ', comment="""{0}"""'.format(self.comment) + string += ')' + return string + + def __reduce__(self): + """ + A helper function used when pickling an ArrheniusEP object. + """ + return (SurfaceArrheniusBM, (self.A, self.n, self.w0, self.E0, self.Tmin, self.Tmax, self.Pmin, self.Pmax, + self.uncertainty, self.comment)) + + cpdef SurfaceArrhenius to_arrhenius(self, double dHrxn): + """ + Return an :class:`Arrhenius` instance of the kinetics model using the + given enthalpy of reaction `dHrxn` to determine the activation energy. + """ + return SurfaceArrhenius( + A=self.A, + n=self.n, + Ea=(self.get_activation_energy(dHrxn) * 0.001, "kJ/mol"), + T0=(1, "K"), + Tmin=self.Tmin, + Tmax=self.Tmax, + uncertainty=self.uncertainty, + comment=self.comment, + ) \ No newline at end of file diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py index dfb598784f..52aa6b4a55 100644 --- a/rmgpy/reaction.py +++ b/rmgpy/reaction.py @@ -53,7 +53,7 @@ from rmgpy.exceptions import ReactionError, KineticsError from rmgpy.kinetics import KineticsData, ArrheniusBM, ArrheniusEP, ThirdBody, Lindemann, Troe, Chebyshev, \ PDepArrhenius, MultiArrhenius, MultiPDepArrhenius, get_rate_coefficient_units_from_reaction_order, \ - SurfaceArrheniusBEP, StickingCoefficientBEP + SurfaceArrheniusBEP, StickingCoefficientBEP, SurfaceArrheniusBM from rmgpy.kinetics.arrhenius import Arrhenius # Separate because we cimport from rmgpy.kinetics.arrhenius from rmgpy.kinetics.surface import SurfaceArrhenius, StickingCoefficient # Separate because we cimport from rmgpy.kinetics.surface from rmgpy.kinetics.diffusionLimited import diffusion_limiter From 58a89954a947496b6b3fdf88ea3aa8ca953f6821 Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 22 Jul 2021 12:04:36 -0400 Subject: [PATCH 05/36] added `add_van_der_waals_bond` mol method --- rmgpy/molecule/molecule.pxd | 2 ++ rmgpy/molecule/molecule.py | 39 +++++++++++++++++++++++++++++++++++++ 2 files changed, 41 insertions(+) diff --git a/rmgpy/molecule/molecule.pxd b/rmgpy/molecule/molecule.pxd index cb720ef84f..12ea1cd9d9 100644 --- a/rmgpy/molecule/molecule.pxd +++ b/rmgpy/molecule/molecule.pxd @@ -170,6 +170,8 @@ cdef class Molecule(Graph): cpdef remove_van_der_waals_bonds(self) + cpdef add_van_der_waals_bond(self) + cpdef sort_atoms(self) cpdef str get_formula(self) diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index d0a057bc2f..b6920296b7 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -1129,6 +1129,45 @@ def remove_van_der_waals_bonds(self): if bond.is_van_der_waals(): self.remove_bond(bond) + def add_van_der_waals_bond(self): + """ + Adds a single van der Waals bond to the graph to connect a surface site to the physisorbed absorbate. + If the adsorbate has lone pairs, the vdw bond will be added to an atom with lone pairs. + Atoms will be selected in the folowing order: + 1) Nitrogen with lone pairs + 2) Oxygen with lone pairs + 3) atoms without lone pairs (C and H) + """ + cython.declare(bond=Bond, fragment=Molecule, adsorbate=Molecule, + surface_site=Molecule,fragments=list, has_surface_site=cython.bint, atom=Atom, lone_pair_atoms=list) + + if self.contains_surface_site(): + fragments = self.split() + if len(fragments) == 2: + has_surface_site = False + for fragment in fragments: + if fragment.is_surface_site(): + has_surface_site = True + surface_site = fragment + else: + adsorbate = fragment + if has_surface_site: + lone_pair_atoms = [atom for atom in adsorbate.atoms if atom.lone_pairs > 0] + if len(lone_pair_atoms) == 1: + # only one atom with a lone pair, so pick this one + atom = lone_pair_atoms[0] + elif len(lone_pair_atoms) > 1: + # more than one atom with lone pair + # sort by element symbol and pick the first (pick nitrogen before oxygen) + lone_pair_atoms.sort(key = lambda atom: atom.symbol) + atom = lone_pair_atoms[0] + else: + # no lone pairs, so pick the first atom + atom = adsorbate.atoms[0] + # create and add the vdw bond + bond = Bond(atom, surface_site.atoms[0],0.0) + self.add_bond(bond) + def sort_atoms(self): """ Sort the atoms in the graph. This can make certain operations, e.g. From eae68570142f05eb7f38e754fe98811cf94205aa Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 22 Jul 2021 12:05:18 -0400 Subject: [PATCH 06/36] added `is_surface_bond` Bond method --- rmgpy/molecule/molecule.pxd | 2 ++ rmgpy/molecule/molecule.py | 10 ++++++++++ 2 files changed, 12 insertions(+) diff --git a/rmgpy/molecule/molecule.pxd b/rmgpy/molecule/molecule.pxd index 12ea1cd9d9..ede13fa00e 100644 --- a/rmgpy/molecule/molecule.pxd +++ b/rmgpy/molecule/molecule.pxd @@ -120,6 +120,8 @@ cdef class Bond(Edge): cpdef bint is_van_der_waals(self) except -2 + cpdef bint is_surface_bond(self) except -2 + cpdef bint is_single(self) except -2 cpdef bint is_double(self) except -2 diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index b6920296b7..772f3b25d9 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -752,6 +752,16 @@ def is_van_der_waals(self): """ return self.is_order(0) + def is_surface_bond(self): + """ + Return ``True`` if the bond represents a surface bond or + ``False`` if not. + """ + if self.atom1.is_surface_site() or self.atom2.is_surface_site(): + return True + else: + return False + def is_order(self, other_order): """ Return ``True`` if the bond is of order other_order or ``False`` if From 512166864736cfbb7d4a9f81377d0998a0a9f21d Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 22 Jul 2021 12:05:47 -0400 Subject: [PATCH 07/36] `get_bde` for metal bonds using MetalDB --- rmgpy/molecule/molecule.py | 67 +++++++++++++++++++++++++++++++++++--- 1 file changed, 63 insertions(+), 4 deletions(-) diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index 772f3b25d9..1e98a15e01 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -45,13 +45,14 @@ import cython import numpy as np +import rmgpy.quantity as quantity import rmgpy.constants as constants import rmgpy.molecule.converter as converter import rmgpy.molecule.element as elements import rmgpy.molecule.group as gr import rmgpy.molecule.resonance as resonance import rmgpy.molecule.translator as translator -from rmgpy.exceptions import DependencyError +from rmgpy.exceptions import DependencyError, DatabaseError from rmgpy.molecule.adjlist import Saturator from rmgpy.molecule.atomtype import AtomType, ATOMTYPES, get_atomtype, AtomTypeError from rmgpy.molecule.element import bdes @@ -638,16 +639,74 @@ def sorting_key(self): self.atom1.number if self.atom1 is not None else 0, self.atom2.number if self.atom2 is not None else 0) - def get_bde(self): + def get_bde(self, metal='Pt111'): """ estimate the bond dissociation energy in J/mol of the bond based on the order of the bond and the atoms involved in the bond + + If the bond involves a surface site, the atomic binding energies from the `metal` (str) + will be use to estimate the bond dissociation energy. """ try: return bdes[(self.atom1.element.symbol, self.atom2.element.symbol, self.order)] except KeyError: - raise KeyError('Bond Dissociation energy not known for combination: ' - '({0},{1},{2})'.format(self.atom1.element.symbol, self.atom2.element.symbol, self.order)) + if not self.is_surface_bond(): + raise KeyError('Bond Dissociation energy not known for combination: ' + '({0},{1},{2})'.format(self.atom1.element.symbol, self.atom2.element.symbol, self.order)) + else: + if self.atom1.is_surface_site(): + adatom = self.atom2.element.symbol + else: + adatom = self.atom1.element.symbol + if self.is_van_der_waals(): + if adatom == 'O': + binding_energy = quantity.Energy(-0.189,'eV/molecule') # binding energy for H2O_ads in surfaceThermoPt111 + elif adatom == 'N': + binding_energy = quantity.Energy(-0.673,'eV/molecule') # binding energy for NH3_ads in surfaceThermoPt111 + elif adatom == 'C': + binding_energy = quantity.Energy(-0.122,'eV/molecule') # binding energy for CH4_ads in surfaceThermoPt111 + elif adatom == 'H': + binding_energy = quantity.Energy(-0.054,'eV/molecule') # binding energy for H2_ads in surfaceThermoPt111 + else: + binding_energy = quantity.Energy(-0.2,'eV/molecule') # default vdw binding energy + return -1 * binding_energy.value_si + try: + from rmgpy.data.rmg import get_db + tdb = get_db('thermo') + metal_db = tdb.surface['metal'] + except (DatabaseError, KeyError): + from rmgpy.data.surface import MetalDatabase + from rmgpy import settings + metal_db = MetalDatabase() + metal_db.load(os.path.join(settings['database.directory'], 'surface')) + metal_binding_energies = metal_db.find_binding_energies(metal=metal) + binding_energy = metal_binding_energies[adatom] + max_bond_order = {'C': 4., 'O': 2., 'N': 3., 'H': 1.} + normalized_bond_order = self.order / max_bond_order[adatom] + # intercepts estimated from DOI:10.1103/PhysRevLett.99.016105 + intercepts = { + 'C' : { + 1 : 0, + 2 : -1.3, + 3 : -1.1, + 4 : 0, + }, + 'N' : { + 1 : -0.5, + 2 : -0.75, + 3 : 0 + }, + 'O' : { + 1 : -0.5, + 2 : 0 + }, + 'H' : { + 1 : 0 + } + } + intercept = quantity.Energy(intercepts[adatom][self.order],'eV/molecule') + binding_energy = quantity.Energy(binding_energy).value_si * normalized_bond_order + intercept.value_si + return -1 * binding_energy def equivalent(self, other): """ From da71141c4b2af2d1d57f9f3c7abfa402745c7e29 Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 22 Jul 2021 12:06:17 -0400 Subject: [PATCH 08/36] added `test_get_bde` unit test --- rmgpy/molecule/moleculeTest.py | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/rmgpy/molecule/moleculeTest.py b/rmgpy/molecule/moleculeTest.py index 576d73c374..0372e51688 100644 --- a/rmgpy/molecule/moleculeTest.py +++ b/rmgpy/molecule/moleculeTest.py @@ -840,6 +840,20 @@ def test_get_bond_string(self): bond = Bond(atom1=Atom(element=get_element(1)), atom2=Atom(element=get_element(6)), order=1) self.assertEqual(bond.get_bond_string(), 'C-H') + def test_get_bde(self): + + bonds = [Bond(atom1=Atom(element=get_element('X')), atom2=Atom(element=get_element(6)), order=0), + Bond(atom1=Atom(element=get_element('X')), atom2=Atom(element=get_element(6)), order=1), + Bond(atom1=Atom(element=get_element('X')), atom2=Atom(element=get_element(6)), order=2), + Bond(atom1=Atom(element=get_element('X')), atom2=Atom(element=get_element(6)), order=3), + Bond(atom1=Atom(element=get_element('X')), atom2=Atom(element=get_element(6)), order=4)] + + bde = 0 + for bond in bonds: + bde_bond = bond.get_bde('Pt111') + self.assertGreater(bde_bond, bde) + bde = bde_bond + ################################################################################ @@ -2759,6 +2773,16 @@ def test_remove_van_der_waals_bonds(self): mol.remove_van_der_waals_bonds() self.assertEqual(len(mol.get_all_edges()), 1) + def test_add_van_der_waals_bond(self): + """Test we can add a van-der-Waals bond""" + mol = Molecule(smiles='*.NOC') + mol.add_van_der_waals_bond() + vdw_bonds = [b for b in mol.get_all_edges() if b.is_van_der_waals()] + self.assertEqual(len(vdw_bonds), 1) + vdw_bond = vdw_bonds[0] + self.assertEqual(vdw_bond.atom1.symbol,'N') + + ################################################################################ if __name__ == '__main__': From d30cad60711ac924f952dd07692f55b4e45fe99b Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 22 Jul 2021 12:06:58 -0400 Subject: [PATCH 09/36] added metal attrs to Reaction --- rmgpy/reaction.pxd | 3 +++ rmgpy/reaction.py | 27 +++++++++++++++++++++++++-- 2 files changed, 28 insertions(+), 2 deletions(-) diff --git a/rmgpy/reaction.pxd b/rmgpy/reaction.pxd index 586c4a6796..77ec42bfc2 100644 --- a/rmgpy/reaction.pxd +++ b/rmgpy/reaction.pxd @@ -58,6 +58,9 @@ cdef class Reaction: cdef public bint is_forward cdef public bint allow_max_rate_violation cdef public object rank + cdef public str metal + cdef public str facet + cdef public str site cpdef bint is_isomerization(self) diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py index 52aa6b4a55..c0345f82d1 100644 --- a/rmgpy/reaction.py +++ b/rmgpy/reaction.py @@ -113,6 +113,9 @@ def __init__(self, allow_max_rate_violation=False, rank=None, comment='', + metal=None, + facet=None, + site=None, is_forward=None, ): self.index = index @@ -134,6 +137,9 @@ def __init__(self, self.is_forward = is_forward self.allow_max_rate_violation = allow_max_rate_violation self.rank = rank + self.metal = metal + self.facet = facet + self.site = site def __repr__(self): """ @@ -157,7 +163,10 @@ def __repr__(self): if self.elementary_high_p: string += 'elementary_high_p={0}, '.format(self.elementary_high_p) if self.comment != '': string += 'comment={0!r}, '.format(self.comment) if self.rank is not None: string += 'rank={0!r},'.format(self.rank) - string = string[:-2] + ')' + if self.metal is not None: string += 'metal={0!r},'.format(self.metal) + if self.facet is not None: string += 'facet={0!r},'.format(self.facet) + if self.site is not None: string += 'site={0!r},'.format(self.site) + string = string[:-1] + ')' return string def __str__(self): @@ -197,8 +206,12 @@ def __reduce__(self): self.pairs, self.allow_pdep_route, self.elementary_high_p, + self.allow_max_rate_violation, self.rank, - self.comment + self.comment, + self.metal, + self.facet, + self.site, )) @property @@ -231,6 +244,16 @@ def degeneracy(self, new): # set new degeneracy self._degeneracy = new + def get_metal_label(self): + """ + retrieve the metal label (str) for the reaction + """ + if self.facet is None: + metal = self.metal # could be None + else: + metal = self.metal + self.facet + return metal + def to_chemkin(self, species_list=None, kinetics=True): """ Return the chemkin-formatted string for this reaction. From 3fcc61ced0018942892f9a2698fe54b40f31bb48 Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 22 Jul 2021 12:07:33 -0400 Subject: [PATCH 10/36] change logging from warning to info --- rmgpy/data/surface.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/rmgpy/data/surface.py b/rmgpy/data/surface.py index a4559adaed..cf52d1662d 100644 --- a/rmgpy/data/surface.py +++ b/rmgpy/data/surface.py @@ -247,9 +247,9 @@ def find_binding_energies(self, metal): metal_binding_energies = self.libraries['surface'].get_binding_energies(metal) except DatabaseError: # no exact match was found, so continue on as if no facet was given - logging.warning("Requested metal %r not found in database.", metal) + logging.info("Requested metal %r not found in database.", metal) metal = metal[:facet.span()[0]] - logging.warning("Searching for generic %r.", metal) + logging.info("Searching for generic %r.", metal) facet = None if facet is None: @@ -263,7 +263,7 @@ def find_binding_energies(self, metal): else: # multiple matches # average the binding energies together? just pick the first one? # just picking the first one for now... - logging.warning(f"Found multiple binding energies for {metal!r}. Using {metal_entry_matches[0]!r}.") + logging.info(f"Found multiple binding energies for {metal!r}. Using {metal_entry_matches[0]!r}.") metal_binding_energies = self.libraries['surface'].get_binding_energies(metal_entry_matches[0]) return metal_binding_energies From 5079f7a495e749641d335cb896989a47facff716 Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 22 Jul 2021 12:08:20 -0400 Subject: [PATCH 11/36] added metal attrs to reaction when loading entries --- rmgpy/data/kinetics/depository.py | 3 ++- rmgpy/data/kinetics/library.py | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/rmgpy/data/kinetics/depository.py b/rmgpy/data/kinetics/depository.py index ad24c7db36..32263ce310 100644 --- a/rmgpy/data/kinetics/depository.py +++ b/rmgpy/data/kinetics/depository.py @@ -221,7 +221,8 @@ def load_entry(self, """ reaction = Reaction(reactants=[], products=[], specific_collider=specificCollider, - degeneracy=degeneracy, duplicate=duplicate, reversible=reversible) + degeneracy=degeneracy, duplicate=duplicate, reversible=reversible, metal=metal, + facet=facet, site=site) entry = Entry( index=index, diff --git a/rmgpy/data/kinetics/library.py b/rmgpy/data/kinetics/library.py index 220b813994..fa479fa731 100644 --- a/rmgpy/data/kinetics/library.py +++ b/rmgpy/data/kinetics/library.py @@ -545,7 +545,7 @@ def load_entry(self, # Make a blank reaction rxn = Reaction(reactants=[], products=[], degeneracy=degeneracy, duplicate=duplicate, reversible=reversible, allow_pdep_route=allow_pdep_route, elementary_high_p=elementary_high_p, - allow_max_rate_violation=allow_max_rate_violation) + allow_max_rate_violation=allow_max_rate_violation,metal=metal, facet=facet, site=site) # if not rxn.is_balanced(): # raise DatabaseError('Reaction {0} in kinetics library {1} was not balanced! Please reformulate.'.format(rxn, self.label)) # label = str(rxn) From f181ca05b8177fb4d26fd37816ed094120954c20 Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 22 Jul 2021 12:09:03 -0400 Subject: [PATCH 12/36] reactions on different metals are not duplicate --- rmgpy/data/kinetics/library.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/rmgpy/data/kinetics/library.py b/rmgpy/data/kinetics/library.py index fa479fa731..fc59c772ac 100644 --- a/rmgpy/data/kinetics/library.py +++ b/rmgpy/data/kinetics/library.py @@ -295,11 +295,16 @@ def check_for_duplicates(self, mark_duplicates=False): """ for entry0 in self.entries.values(): reaction0 = entry0.item + metal0 = entry0.get_metal_label() if not reaction0.duplicate: # This reaction is not marked as a duplicate reaction # This means that if we find any duplicate reactions, it is an error for entry in self.entries.values(): reaction = entry.item + metal = entry.get_metal_label() + if metal0 != metal: + # the metals don't match, so they are not duplicate reactions + continue if reaction0 is not reaction and reaction0.is_isomorphic(reaction): # We found a duplicate reaction that wasn't marked! # RMG requires all duplicate reactions to be marked, unlike CHEMKIN From 4a66dc469cf5416f811251de2995d8a2b4ae466e Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 22 Jul 2021 12:16:18 -0400 Subject: [PATCH 13/36] added metal and facet to input, RMG and CoreEdgeReactionModel classes this is needed so we can select reactions that match metals in training and reaction libraries --- rmgpy/rmg/input.py | 10 +++++++--- rmgpy/rmg/main.py | 9 ++++++++- rmgpy/rmg/model.py | 2 ++ 3 files changed, 17 insertions(+), 4 deletions(-) diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py index 181024457e..3c4fb215ec 100644 --- a/rmgpy/rmg/input.py +++ b/rmgpy/rmg/input.py @@ -101,6 +101,7 @@ def database( def catalyst_properties(bindingEnergies=None, surfaceSiteDensity=None, metal=None, + facet=None, coverageDependence=False): """ Specify the properties of the catalyst. @@ -121,10 +122,13 @@ def catalyst_properties(bindingEnergies=None, if metal: try: logging.info("Using catalyst surface properties from metal %r.", metal) - rmg.binding_energies = metal_db.get_binding_energies(metal) - rmg.surface_site_density = metal_db.get_surface_site_density(metal) + rmg.metal = metal + rmg.facet = facet + + rmg.binding_energies = metal_db.get_binding_energies(metal+facet) + rmg.surface_site_density = metal_db.get_surface_site_density(metal+facet) except DatabaseError: - logging.error('Metal %r missing from surface library. Please specify both metal and facet.', metal) + logging.error('Metal %r missing from surface library. Please specify both metal and facet.', metal+facet) raise else: # metal not specified if bindingEnergies is None: diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index bd7f8f896d..aa847d4846 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -184,6 +184,8 @@ def clear(self): self.kinetics_estimator = 'group additivity' self.solvent = None self.diffusion_limiter = None + self.metal = None + self.facet = None self.surface_site_density = None self.binding_energies = None self.coverage_dependence = False @@ -268,7 +270,12 @@ def load_input(self, path=None): self.reaction_model.core.phase_system.phases["Surface"].site_density = self.surface_site_density.value_si self.reaction_model.edge.phase_system.phases["Surface"].site_density = self.surface_site_density.value_si self.reaction_model.coverage_dependence = self.coverage_dependence - + + if self.metal: + self.reaction_model.metal = self.metal + if self.facet: + self.reaction_model.facet = self.facet + self.reaction_model.verbose_comments = self.verbose_comments self.reaction_model.save_edge_species = self.save_edge_species diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index 90d73ea989..a3bc0c319b 100644 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -235,6 +235,8 @@ def __init__(self, core=None, edge=None, surface=None): 3 R!H u1 px c0 {2,S} 4 H u0 p0 c0 {1,S} """)] + self.metal = None + self.facet = None def check_for_existing_species(self, molecule): """ From a3fd3e9d96834b70216f62268798e67dcdaf9335 Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 22 Jul 2021 12:19:13 -0400 Subject: [PATCH 14/36] select the best kinetics with metal and facet attrs --- rmgpy/data/kinetics/family.py | 17 +++++++++++++++-- rmgpy/rmg/model.py | 3 +++ 2 files changed, 18 insertions(+), 2 deletions(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 0a7b14ca2e..f04c7e8234 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -2660,12 +2660,22 @@ def get_kinetics_from_depository(self, depository, reaction, template, degenerac kinetics.comment += "\nsite: {}".format(self.site) return kinetics_list - def _select_best_kinetics(self, kinetics_list): + def _select_best_kinetics(self, kinetics_list, metal=None, facet=None): """ For a given set of kinetics `kinetics_list`, return the kinetics deemed to be the "best". This is determined to be the one with the lowest non-zero rank that occurs first (has the lowest index). """ + if metal: + # select only the reactions that match metals + _kinetics_list = [k for k in kinetics_list if k[1].metal == metal] + if _kinetics_list: + kinetics_list = _kinetics_list + if facet: + # select only the reactions that match facets + _kinetics_list = [k for k in kinetics_list if k[1].facet == facet] + if _kinetics_list: + kinetics_list = _kinetics_list if any([x[1].rank == 0 for x in kinetics_list]) and not all([x[1].rank == 0 for x in kinetics_list]): kinetics_list = [x for x in kinetics_list if x[1].rank != 0] kinetics_list.sort(key=lambda x: (x[1].rank, x[1].index)) @@ -2695,11 +2705,14 @@ def get_kinetics(self, reaction, template_labels, degeneracy=1, estimator='', re template = self.retrieve_template(template_labels) + metal = reaction.metal + facet = reaction.facet + # Check the various depositories for kinetics for depository in depositories: kinetics_list0 = self.get_kinetics_from_depository(depository, reaction, template, degeneracy) if len(kinetics_list0) > 0 and not return_all_kinetics: - kinetics, entry, is_forward = self._select_best_kinetics(kinetics_list0) + kinetics, entry, is_forward = self._select_best_kinetics(kinetics_list0, metal=metal, facet=facet) return kinetics, depository, entry, is_forward else: for kinetics, entry, is_forward in kinetics_list0: diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index a3bc0c319b..a5c479035a 100644 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -920,6 +920,9 @@ def generate_kinetics(self, reaction): assert isinstance(reaction, TemplateReaction) family = get_family_library_object(reaction.family) + if reaction.is_surface_reaction(): + reaction.metal = self.metal + reaction.facet = self.facet # Get the kinetics for the reaction kinetics, source, entry, is_forward = family.get_kinetics(reaction, template_labels=reaction.template, From 3d2daa815e5f75069309c0949ed7dd4c73989a88 Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 22 Jul 2021 12:24:26 -0400 Subject: [PATCH 15/36] added `vdw_bonds` family attr this is needed for autogen trees with vdw bonds in the template in order to avoid splitting the vdw bonded complex into difference reactants --- rmgpy/data/kinetics/family.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index f04c7e8234..518de92d4a 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -727,6 +727,7 @@ def load(self, path, local_context=None, global_context=None, depository_labels= self.boundary_atoms = local_context.get('boundaryAtoms', None) self.tree_distances = local_context.get('treeDistances', None) self.reverse_map = local_context.get('reverseMap', None) + self.vdw_bonds = local_context.get('vdwBonds', None) self.reactant_num = local_context.get('reactantNum', None) self.product_num = local_context.get('productNum', None) @@ -1014,6 +1015,9 @@ def save_groups(self, path): if self.reverse_map is not None: f.write('reverseMap = {0}\n\n'.format(self.reverse_map)) + if self.vdw_bonds is not None: + f.write('vdwBonds = {0}\n\n'.format(self.vdw_bonds)) + if self.reactant_num is not None: f.write('reactantNum = {0}\n\n'.format(self.reactant_num)) if self.product_num is not None: From c209e60fe7484932b633b1e8fbb72f5681e053fc Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 22 Jul 2021 12:33:02 -0400 Subject: [PATCH 16/36] added `split_template` family method this method is needed to split merged templates (generated in autotrees) for multi-molecular families when generating reactions --- rmgpy/data/kinetics/family.py | 63 ++++++++++++++++++++--------------- rmgpy/rmg/main.py | 1 + 2 files changed, 37 insertions(+), 27 deletions(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 518de92d4a..b0312785bb 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -859,6 +859,35 @@ def save_entry(self, f, entry): """ return save_entry(f, entry) + def split_template(self): + """ + If the family has one template and is not unimolecular, + split the reaction template into multiple reactants + """ + if self.reactant_num > 1 and len(self.forward_template.reactants) == 1: + template = self.forward_template.reactants[0] + if isinstance(template.item, Group): + if self.vdw_bonds: + for label0,label1 in self.vdw_bonds.items(): + atom0 = template.item.get_labeled_atoms(label0)[0] + atom1 = template.item.get_labeled_atoms(label1)[0] + bond = GroupBond(atom0, atom1, order=[0]) + template.item.add_bond(bond) + + groups = template.item.split() + template_reactants = [] + if len(groups) > 1: + for i,group in enumerate(groups): + group.remove_van_der_waals_bonds() + template_reactant = deepcopy(template) + template_reactant.item = group + template_reactant.label += str(i) + template_reactants.append(template_reactant) + self.forward_template.reactants = template_reactants + if self.reverse_template: + self.reverse_template.products = template_reactants + assert len(self.forward_template.reactants) == self.reactant_num + def save_training_reactions(self, reactions, reference=None, reference_type='', short_desc='', long_desc='', rank=3): """ @@ -2064,6 +2093,7 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson # This also makes a copy of the reactants list so we don't modify the # original reactants = [reactant if isinstance(reactant, list) else [reactant] for reactant in reactants] + self.split_template() if forward: template = self.forward_template @@ -2077,21 +2107,7 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson if self.auto_generated and reactant_num != len(reactants): return [] - if len(reactants) > len(template.reactants): - # If the template contains a surface site, we do not want to split it because it will break vdw bonds - if isinstance(template.reactants[0].item, Group): - if template.reactants[0].item.contains_surface_site(): - return [] - # if the family has one template and is bimolecular split template into multiple reactants - try: - grps = template.reactants[0].item.split() - template_reactants = [] - for grp in grps: - template_reactants.append(grp) - except AttributeError: - template_reactants = [x.item for x in template.reactants] - else: - template_reactants = [x.item for x in template.reactants] + template_reactants = [x.item for x in template.reactants] # Unimolecular reactants: A --> products if len(reactants) == 1 and len(template_reactants) == 1: @@ -2831,23 +2847,16 @@ def get_labeled_reactants_and_products(self, reactants, products): new entities in memory so input molecules `reactants` and `products` won't be affected. If RMG cannot find appropriate labels, (None, None) will be returned. """ + + self.split_template() # if the family has one template and is bimolecular split template into multiple reactants template = self.forward_template reactants0 = [reactant.copy(deep=True) for reactant in reactants] + template_reactants = [x.item for x in template.reactants] if self.auto_generated and self.reactant_num != len(reactants): return None, None - - if len(reactants) > len(template.reactants): - # if the family has one template and is bimolecular split template into multiple reactants - try: - grps = template.reactants[0].item.split() - template_reactants = [] - for grp in grps: - template_reactants.append(grp) - except AttributeError: - template_reactants = [x.item for x in template.reactants] - else: - template_reactants = [x.item for x in template.reactants] + + template_reactants = [x.item for x in template.reactants] if len(reactants0) == 1: molecule = reactants0[0] diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index aa847d4846..c4676fd934 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -455,6 +455,7 @@ def load_database(self): for family in self.database.kinetics.families.values(): if not family.auto_generated: family.fill_rules_by_averaging_up(verbose=self.verbose_comments) + family.split_template() def initialize(self, **kwargs): """ From 9fdc4be346d70fae12efb8649e11f1c9526852ec Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 22 Jul 2021 12:34:49 -0400 Subject: [PATCH 17/36] use SurfaceArrheniusBM for surace reactions with autotree gen methods --- rmgpy/data/kinetics/family.py | 53 +++++++++++++++++++++++++---------- 1 file changed, 38 insertions(+), 15 deletions(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index b0312785bb..a1da40599d 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -3734,6 +3734,11 @@ def cross_validate(self, folds=5, template_rxn_map=None, test_rxn_inds=None, T=1 if template_rxn_map is None: template_rxn_map = self.get_reaction_matches(remove_degeneracy=True, get_reverse=True, fix_labels=True) + if template_rxn_map['Root'][0].is_surface_reaction(): + surface = True + else: + surface = False + rxns = np.array(template_rxn_map['Root']) if test_rxn_inds is None: @@ -3785,7 +3790,10 @@ def cross_validate(self, folds=5, template_rxn_map=None, test_rxn_inds=None, T=1 L = list(set(template_rxn_map[entry.label]) - set(rxns_test)) if L != []: - kinetics = ArrheniusBM().fit_to_reactions(L, recipe=self.forward_recipe.actions) + if surface: + kinetics = SurfaceArrheniusBM().fit_to_reactions(L, recipe=self.forward_recipe.actions) + else: + kinetics = ArrheniusBM().fit_to_reactions(L, recipe=self.forward_recipe.actions) kinetics = kinetics.to_arrhenius(rxn.get_enthalpy_of_reaction(T)) k = kinetics.get_rate_coefficient(T) errors[rxn] = np.log(k / krxn) @@ -4689,14 +4697,31 @@ def _make_rule(rr): """ recipe, rxns, Tref, fmax, label, ranks = rr n = len(rxns) + if n == 0: + return None for i, rxn in enumerate(rxns): rxn.rank = ranks[i] + if rxn.is_surface_reaction(): + surface = True + else: + surface = False rxns = np.array(rxns) data_mean = np.mean(np.log([r.kinetics.get_rate_coefficient(Tref) for r in rxns])) - if n > 0: + if surface: + kin = SurfaceArrheniusBM().fit_to_reactions(rxns, recipe=recipe) + else: kin = ArrheniusBM().fit_to_reactions(rxns, recipe=recipe) - if n == 1: - kin.uncertainty = RateUncertainty(mu=0.0, var=(np.log(fmax) / 2.0) ** 2, N=1, Tref=Tref, data_mean=data_mean, correlation=label) + if n == 1: + kin.uncertainty = RateUncertainty(mu=0.0, var=(np.log(fmax) / 2.0) ** 2, N=1, Tref=Tref, data_mean=data_mean, correlation=label) + else: + if surface: + dlnks = np.array([ + np.log( + SurfaceArrheniusBM().fit_to_reactions(rxns[list(set(range(len(rxns))) - {i})], recipe=recipe) + .to_arrhenius(rxn.get_enthalpy_of_reaction(Tref)) + .get_rate_coefficient(T=Tref) / rxn.get_rate_coefficient(T=Tref) + ) for i, rxn in enumerate(rxns) + ]) # 1) fit to set of reactions without the current reaction (k) 2) compute log(kfit/kactual) at Tref else: dlnks = np.array([ np.log( @@ -4705,17 +4730,15 @@ def _make_rule(rr): .get_rate_coefficient(T=Tref) / rxn.get_rate_coefficient(T=Tref) ) for i, rxn in enumerate(rxns) ]) # 1) fit to set of reactions without the current reaction (k) 2) compute log(kfit/kactual) at Tref - varis = (np.array([rank_accuracy_map[rxn.rank].value_si for rxn in rxns]) / (2.0 * 8.314 * Tref)) ** 2 - # weighted average calculations - ws = 1.0 / varis - V1 = ws.sum() - V2 = (ws ** 2).sum() - mu = np.dot(ws, dlnks) / V1 - s = np.sqrt(np.dot(ws, (dlnks - mu) ** 2) / (V1 - V2 / V1)) - kin.uncertainty = RateUncertainty(mu=mu, var=s ** 2, N=n, Tref=Tref, data_mean=data_mean, correlation=label) - return kin - else: - return None + varis = (np.array([rank_accuracy_map[rxn.rank].value_si for rxn in rxns]) / (2.0 * 8.314 * Tref)) ** 2 + # weighted average calculations + ws = 1.0 / varis + V1 = ws.sum() + V2 = (ws ** 2).sum() + mu = np.dot(ws, dlnks) / V1 + s = np.sqrt(np.dot(ws, (dlnks - mu) ** 2) / (V1 - V2 / V1)) + kin.uncertainty = RateUncertainty(mu=mu, var=s ** 2, N=n, Tref=Tref, data_mean=data_mean, correlation=label) + return kin def _spawn_tree_process(family, template_rxn_map, obj, T, nprocs, depth, min_splitable_entry_num, min_rxns_to_spawn, extension_iter_max, extension_iter_item_cap): From f75f7c8ac9f3b50653b40462b4df19d5904a0184 Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 22 Jul 2021 12:35:50 -0400 Subject: [PATCH 18/36] save training reactions with metal attrs --- rmgpy/data/kinetics/family.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index a1da40599d..37fc22bab1 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -889,7 +889,7 @@ def split_template(self): assert len(self.forward_template.reactants) == self.reactant_num def save_training_reactions(self, reactions, reference=None, reference_type='', short_desc='', long_desc='', - rank=3): + metal=None, facet=None, site=None, rank=3): """ This function takes a list of reactions appends it to the training reactions file. It ignores the existence of duplicate reactions. @@ -908,6 +908,12 @@ def save_training_reactions(self, reactions, reference=None, reference_type='', short_desc = [short_desc] * len(reactions) if not isinstance(long_desc, list): long_desc = [long_desc] * len(reactions) + if not isinstance(metal, list): + metal = [metal] * len(reactions) + if not isinstance(facet, list): + facet = [facet] * len(reactions) + if not isinstance(site, list): + site = [site] * len(reactions) if not isinstance(rank, list): rank = [rank] * len(reactions) @@ -978,6 +984,9 @@ def save_training_reactions(self, reactions, reference=None, reference_type='', short_desc=str(short_desc[i]), long_desc=str(long_desc[i]), rank=rank[i], + metal=metal[i], + facet=facet[i], + site=site[i] ) # Add this entry to the loaded depository so it is immediately usable From dbab14b7027348c0ffdbac96aad0b4c6a2b3267a Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 22 Jul 2021 12:37:34 -0400 Subject: [PATCH 19/36] set family reactant_num to length of template if not specified --- rmgpy/data/kinetics/family.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 37fc22bab1..cbce87b31e 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -734,10 +734,9 @@ def load(self, path, local_context=None, global_context=None, depository_labels= self.auto_generated = local_context.get('autoGenerated', False) - if self.reactant_num: - self.groups.reactant_num = self.reactant_num - else: - self.groups.reactant_num = len(self.forward_template.reactants) + if not self.reactant_num: + self.reactant_num = len(self.forward_template.reactants) + self.groups.reactant_num = self.reactant_num # Generate the reverse template if necessary self.forward_template.reactants = [self.groups.entries[label] for label in self.forward_template.reactants] From cae9258ef7193a350d7a4676af65f1d05b2ff5d2 Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 22 Jul 2021 12:38:40 -0400 Subject: [PATCH 20/36] using get_metal_label function --- rmgpy/data/kinetics/family.py | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index cbce87b31e..b8d6ed7bc6 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -1313,10 +1313,7 @@ def add_rules_from_training(self, thermo_database=None, train_indices=None): if quantum_mechanics: quantum_mechanics.run_jobs(item.reactants + item.products, procnum=procnum) - if entry.facet is None: - metal = entry.metal # could be None - else: - metal = entry.metal + entry.facet + metal = entry.get_metal_label() for reactant in item.reactants: # Clear atom labels to avoid effects on thermo generation, ok because this is a deepcopy @@ -4204,10 +4201,7 @@ def get_reactant_thermo(reactant,metal): for i, entry in enumerate(entries): if estimate_thermo: # parse out the metal to scale to - if entry.facet is None: - metal = entry.metal # could be None - else: - metal = entry.metal + entry.facet + metal = entry.get_metal_label() for j, react in enumerate(entry.item.reactants): if rxns[i].reactants[j].thermo is None: rxns[i].reactants[j].thermo = get_reactant_thermo(react,metal) From ab9d98148e47453efdbfabc27c9507c8f6383abe Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 22 Jul 2021 12:41:12 -0400 Subject: [PATCH 21/36] get_training_set for surface families - make deepcopies of adsorbates since thermo is different on different metals - assign metal/facet/site to reverse entries - stick entry index on reaction index (helpful for debugging) --- rmgpy/data/kinetics/family.py | 37 ++++++++++++++++++++++++++--------- 1 file changed, 28 insertions(+), 9 deletions(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index b8d6ed7bc6..6147579cb6 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -4203,14 +4203,23 @@ def get_reactant_thermo(reactant,metal): # parse out the metal to scale to metal = entry.get_metal_label() for j, react in enumerate(entry.item.reactants): + if metal: + # we need to make a new species for different metals beacause thermo is different + rxns[i].reactants[j] = react.copy(deep=True) + rxns[i].reactants[j].thermo = None if rxns[i].reactants[j].thermo is None: - rxns[i].reactants[j].thermo = get_reactant_thermo(react,metal) + rxns[i].reactants[j].thermo = get_reactant_thermo(rxns[i].reactants[j],metal) for j, react in enumerate(entry.item.products): + if metal: + # we need to make a new species for different metals beacause thermo is different + rxns[i].products[j] = react.copy(deep=True) + rxns[i].products[j].thermo = None if rxns[i].products[j].thermo is None: - rxns[i].products[j].thermo = get_reactant_thermo(react,metal) + rxns[i].products[j].thermo = get_reactant_thermo(rxns[i].products[j],metal) rxns[i].kinetics = entry.data rxns[i].rank = entry.rank + rxns[i].index = entry.index if remove_degeneracy: # adjust for degeneracy rxns[i].kinetics.A.value_si /= rxns[i].degeneracy @@ -4285,9 +4294,13 @@ def get_reactant_thermo(reactant,metal): rrev.is_forward = False if estimate_thermo: - for rev_react in rrev.reactants: - if rev_react.thermo is None: - rev_react.thermo = get_reactant_thermo(rev_react,metal) + for x,react in enumerate(rrev.reactants): + if metal: + # we need to make a new species for different metals beacause thermo is different + rrev.reactants[x] = react.copy(deep=True) + rrev.reactants[x].thermo = None + if rrev.reactants[x].thermo is None: + rrev.reactants[x].thermo = get_reactant_thermo(rrev.reactants[x],metal) rev_rxns.append(rrev) @@ -4319,15 +4332,21 @@ def get_reactant_thermo(reactant,metal): products = [Species(molecule=[p]) for p in products] rrev = Reaction(reactants=products, products=rxns[i].reactants, - kinetics=rxns[i].generate_reverse_rate_coefficient(), rank=rxns[i].rank) + kinetics=rxns[i].generate_reverse_rate_coefficient(), rank=rxns[i].rank, + metal = rxns[i].metal, facet = rxns[i].facet, site = rxns[i].site) rrev.is_forward = False if estimate_thermo: - for rev_react in rrev.reactants: - if rev_react.thermo is None: - rev_react.thermo = get_reactant_thermo(rev_react,metal) + for x,react in enumerate(rrev.reactants): + if metal: + # we need to make a new species for different metals beacause thermo is different + rrev.reactants[x] = react.copy(deep=True) + rrev.reactants[x].thermo = None + if rrev.reactants[x].thermo is None: + rrev.reactants[x].thermo = get_reactant_thermo(rrev.reactants[x],metal) rxns[i] = rrev + rxns[i].index = entry.index if self.own_reverse and get_reverse: return rxns + rev_rxns From 468331ef5c0ac437e0536deeb62e6067dd37fc68 Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 22 Jul 2021 12:41:39 -0400 Subject: [PATCH 22/36] `get_kinetics_from_depository` fixup --- rmgpy/data/kinetics/family.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 6147579cb6..d0b47b712d 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -2678,11 +2678,11 @@ def get_kinetics_from_depository(self, depository, reaction, template, degenerac '[{0}]'.format(';'.join([g.label for g in template]))) kinetics.comment += "\nfamily: {}".format(self.label) if entry.metal: - kinetics.comment += "\nmetal: {}".format(self.metal) + kinetics.comment += "\nmetal: {}".format(entry.metal) if entry.facet: - kinetics.comment += "\nfacet: {}".format(self.facet) - if entry.facet: - kinetics.comment += "\nsite: {}".format(self.site) + kinetics.comment += "\nfacet: {}".format(entry.facet) + if entry.site: + kinetics.comment += "\nsite: {}".format(entry.site) return kinetics_list def _select_best_kinetics(self, kinetics_list, metal=None, facet=None): From 675c18ae9ceef2d58081509099e8594581e2fa48 Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 22 Jul 2021 17:39:46 -0400 Subject: [PATCH 23/36] added metal attrs to kinetics database and generate reactions metal attrs are loaded from input and set as attrs on kinetics database. They are taken from kinetics database and set on reactions generated with the kineticsDB --- rmgpy/data/kinetics/database.py | 7 ++++++- rmgpy/data/kinetics/family.py | 31 ++++++++++++++++++++++++++----- rmgpy/rmg/input.py | 12 ++++++++---- rmgpy/rmg/main.py | 4 ++++ rmgpy/rmg/model.py | 3 --- 5 files changed, 44 insertions(+), 13 deletions(-) diff --git a/rmgpy/data/kinetics/database.py b/rmgpy/data/kinetics/database.py index 7525f5ed62..aebb83568a 100644 --- a/rmgpy/data/kinetics/database.py +++ b/rmgpy/data/kinetics/database.py @@ -84,6 +84,9 @@ def __init__(self): 'SurfaceArrheniusBM': SurfaceArrheniusBM } self.global_context = {} + self.metal = None + self.facet = None + self.site = None def __reduce__(self): """ @@ -560,7 +563,9 @@ def react_molecules(self, molecules, products=None, only_families=None, prod_res if only_families is None or label in only_families: try: reaction_list.extend(family.generate_reactions(molecules, products=products, - prod_resonance=prod_resonance)) + prod_resonance=prod_resonance, + metal=self.metal, facet=self.facet, + site=self.site)) except: logging.error("Problem family: {}".format(label)) logging.error("Problem reactants: {}".format(molecules)) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index d0b47b712d..669b6f5f6f 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -1850,7 +1850,8 @@ def _match_reactant_to_template(self, reactant, template_reactant): else: raise NotImplementedError("Not expecting template of type {}".format(type(struct))) - def generate_reactions(self, reactants, products=None, prod_resonance=True, delete_labels=True, relabel_atoms=True): + def generate_reactions(self, reactants, products=None, prod_resonance=True, delete_labels=True, relabel_atoms=True, + metal= None, facet=None, site=None): """ Generate all reactions between the provided list of one, two, or three `reactants`, which should be either single :class:`Molecule` objects @@ -1885,6 +1886,9 @@ def generate_reactions(self, reactants, products=None, prod_resonance=True, dele prod_resonance=prod_resonance, delete_labels=delete_labels, relabel_atoms=relabel_atoms, + metal=metal, + facet=facet, + site=site )) if not self.own_reverse and self.reversible: @@ -1896,6 +1900,9 @@ def generate_reactions(self, reactants, products=None, prod_resonance=True, dele prod_resonance=prod_resonance, delete_labels=delete_labels, relabel_atoms=relabel_atoms, + metal=metal, + facet=facet, + site=site )) return reaction_list @@ -1907,6 +1914,9 @@ def add_reverse_attribute(self, rxn, react_non_reactive=True): Returns `True` if successful and `False` if the reverse reaction is forbidden. Will raise a `KineticsError` if unsuccessful for other reasons. """ + + #parse out metal attrs of rxn + metal, facet, site = rxn.metal, rxn.facet, rxn.site if self.own_reverse and all([spc.has_reactive_molecule() for spc in rxn.products]): # Check if the reactants are the same same_reactants = 0 @@ -1926,7 +1936,8 @@ def add_reverse_attribute(self, rxn, react_non_reactive=True): reaction_list = self._generate_reactions([spc.molecule for spc in rxn.products], products=rxn.reactants, forward=True, - react_non_reactive=react_non_reactive) + react_non_reactive=react_non_reactive, + metal=metal, facet=facet, site=site) reactions = find_degenerate_reactions(reaction_list, same_reactants, kinetics_family=self) if len(reactions) == 0: logging.error("Expecting one matching reverse reaction, not zero in reaction family {0} for " @@ -1949,7 +1960,8 @@ def add_reverse_attribute(self, rxn, react_non_reactive=True): try: reaction_list = self._generate_reactions([spc.molecule for spc in rxn.products], products=rxn.reactants, forward=True, - react_non_reactive=react_non_reactive) + react_non_reactive=react_non_reactive, + metal=metal, facet=facet, site=site) reactions = find_degenerate_reactions(reaction_list, same_reactants, kinetics_family=self) finally: self.forbidden = temp_object @@ -1998,6 +2010,8 @@ def calculate_degeneracy(self, reaction): identical reactants (since you will be reducing them later in the algorithm), add `ignoreSameReactants= True` to this method. """ + # parse out metal attrs of reaction + metal, facet, site = reaction.metal, reaction.facet, reaction.site # 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 @@ -2042,7 +2056,7 @@ def calculate_degeneracy(self, reaction): reactions = [] for combo in molecule_combos: reactions.extend(self._generate_reactions(combo, products=reaction.products, forward=True, - react_non_reactive=True)) + react_non_reactive=True, metal=metal, facet=facet, site=site)) # remove degenerate reactions reactions = find_degenerate_reactions(reactions, same_reactants, template=reaction.template, @@ -2060,7 +2074,8 @@ def calculate_degeneracy(self, reaction): return reactions[0].degeneracy def _generate_reactions(self, reactants, products=None, forward=True, prod_resonance=True, - react_non_reactive=False, delete_labels=True, relabel_atoms=True): + react_non_reactive=False, delete_labels=True, relabel_atoms=True, + metal=None, facet=None, site=None): """ Generate a list of all the possible reactions of this family between the list of `reactants`. The number of reactants provided must match @@ -2434,6 +2449,12 @@ def generate_products_and_reactions(order): # Also store the reaction template (useful so we can easily get the kinetics later) for reaction in rxn_list: + if any([metal,facet,site]): + if reaction.is_surface_reaction(): + reaction.metal = metal + reaction.facet = facet + reaction.site = site + # Restore the labeled atoms long enough to generate some metadata for reactant in reaction.reactants: reactant.clear_labeled_atoms() diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py index 3c4fb215ec..8521a9b596 100644 --- a/rmgpy/rmg/input.py +++ b/rmgpy/rmg/input.py @@ -120,15 +120,19 @@ def catalyst_properties(bindingEnergies=None, "or the surfaceSiteDensity and bindingEnergies, but not both.") if metal: + if facet: + metal_label = metal + facet + else: + metal_label = metal try: - logging.info("Using catalyst surface properties from metal %r.", metal) + logging.info("Using catalyst surface properties from metal %r.", metal_label) rmg.metal = metal rmg.facet = facet - rmg.binding_energies = metal_db.get_binding_energies(metal+facet) - rmg.surface_site_density = metal_db.get_surface_site_density(metal+facet) + rmg.binding_energies = metal_db.get_binding_energies(metal_label) + rmg.surface_site_density = metal_db.get_surface_site_density(metal_label) except DatabaseError: - logging.error('Metal %r missing from surface library. Please specify both metal and facet.', metal+facet) + logging.error('Metal %r missing from surface library. Please specify both metal and facet.', metal_label) raise else: # metal not specified if bindingEnergies is None: diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index c4676fd934..8e5a706a49 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -389,6 +389,10 @@ def load_database(self): depository=False, # Don't bother loading the depository information, as we don't use it ) + # add metal attrs to kinetics database + self.database.kinetics.metal = self.metal + self.database.kinetics.facet = self.facet + # Turn off reversibility for families with three products if desired if not self.trimolecular_product_reversible: for family in self.database.kinetics.families.values(): diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index a5c479035a..a3bc0c319b 100644 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -920,9 +920,6 @@ def generate_kinetics(self, reaction): assert isinstance(reaction, TemplateReaction) family = get_family_library_object(reaction.family) - if reaction.is_surface_reaction(): - reaction.metal = self.metal - reaction.facet = self.facet # Get the kinetics for the reaction kinetics, source, entry, is_forward = family.get_kinetics(reaction, template_labels=reaction.template, From e4822a914aea53abdcaf3b79efe2c0ae253ecc6e Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 22 Jul 2021 18:10:06 -0400 Subject: [PATCH 24/36] added metal attrs to kinetics Groups and descend tree algorithm --- rmgpy/data/base.py | 63 ++++++++++++++++++++++++++--------- rmgpy/data/baseTest.py | 33 ++++++++++++++++++ rmgpy/data/kinetics/family.py | 25 ++++++++++---- rmgpy/data/kinetics/groups.py | 12 ++++--- 4 files changed, 107 insertions(+), 26 deletions(-) diff --git a/rmgpy/data/base.py b/rmgpy/data/base.py index 5850bd3bad..94257f7369 100644 --- a/rmgpy/data/base.py +++ b/rmgpy/data/base.py @@ -911,8 +911,10 @@ def match_node_to_node(self, node, node_other): Both `node` and `node_other` must be Entry types with items containing Group or LogicNode types. """ if isinstance(node.item, Group) and isinstance(node_other.item, Group): - return (self.match_node_to_structure(node, node_other.item, atoms=node_other.item.get_all_labeled_atoms()) and - self.match_node_to_structure(node_other, node.item, atoms=node.item.get_all_labeled_atoms())) + return (self.match_node_to_structure(node, node_other.item, atoms=node_other.item.get_all_labeled_atoms(), + metal=node_other.metal, facet=node_other.facet, site=node_other.site) and + self.match_node_to_structure(node_other, node.item, atoms=node.item.get_all_labeled_atoms(), + metal=node.metal, facet=node.facet, site=node.site)) elif isinstance(node.item, LogicOr) and isinstance(node_other.item, LogicOr): return node.item.match_logic_or(node_other.item) else: @@ -928,9 +930,11 @@ def match_node_to_child(self, parent_node, child_node): if isinstance(parent_node.item, Group) and isinstance(child_node.item, Group): if (self.match_node_to_structure(parent_node, child_node.item, - atoms=child_node.item.get_all_labeled_atoms(), strict=True) and + atoms=child_node.item.get_all_labeled_atoms(), strict=True, + metal=child_node.metal, facet=child_node.facet, site=child_node.site) and not self.match_node_to_structure(child_node, parent_node.item, - atoms=parent_node.item.get_all_labeled_atoms(), strict=True)): + atoms=parent_node.item.get_all_labeled_atoms(), strict=True, + metal=parent_node.metal, facet=parent_node.facet, site=parent_node.site)): return True return False @@ -950,7 +954,7 @@ def match_node_to_child(self, parent_node, child_node): elif isinstance(parent_node.item, LogicOr): return child_node.label in parent_node.item.components - def match_node_to_structure(self, node, structure, atoms, strict=False): + def match_node_to_structure(self, node, structure, atoms, strict=False, metal=None, facet=None, site=None): """ Return :data:`True` if the `structure` centered at `atom` matches the structure at `node` in the dictionary. The structure at `node` should @@ -963,6 +967,9 @@ def match_node_to_structure(self, node, structure, atoms, strict=False): Matching to structure is more strict than to node. All labels in structure must be found in node. However the reverse is not true, unless `strict` is set to True. + + If the node has a metal attribute (metal, facet, and/or site) that does not match the + provdied `metal`, `facet`, and/or `site`, the structure is not a match (returns `False`) =================== ======================================================== Attribute Description @@ -971,11 +978,25 @@ def match_node_to_structure(self, node, structure, atoms, strict=False): `structure` A Group or a Molecule `atoms` Dictionary of {label: atom} in the structure. A possible dictionary is the one produced by structure.get_all_labeled_atoms() `strict` If set to ``True``, ensures that all the node's atomLabels are matched by in the structure + `metal` metal of structure (default is None) + `facet` facet of structure (default is None) + `site` site of structure (default is None) =================== ======================================================== """ if isinstance(node, str): node = self.entries[node] group = node.item + + if node.metal: + if metal != node.metal: + return False + if node.facet: + if facet != node.facet: + return False + if node.site: + if site != node.site: + return False + if isinstance(group, LogicNode): return group.match_to_structure(self, structure, atoms, strict) else: @@ -1033,7 +1054,7 @@ def match_node_to_structure(self, node, structure, atoms, strict=False): return result - def descend_tree(self, structure, atoms, root=None, strict=False): + def descend_tree(self, structure, atoms, root=None, strict=False, metal=None, facet=None, site=None): """ Descend the tree in search of the functional group node that best matches the local structure around `atoms` in `structure`. @@ -1045,24 +1066,28 @@ def descend_tree(self, structure, atoms, root=None, strict=False): Set strict to ``True`` if all labels in final matched node must match that of the structure. This is used in kinetics groups to find the correct reaction template, but not generally used in other GAVs due to species generally not being prelabeled. + + `metal` (None or str): metal of structure (default is None) + `facet` (None or str): facet of structure (default is None) + `site` (None or str): site of structure (default is None) """ if root is None: for root in self.top: - if self.match_node_to_structure(root, structure, atoms, strict): + if self.match_node_to_structure(root, structure, atoms, strict, metal=metal, facet=facet, site=site): break # We've found a matching root else: # didn't break - matched no top nodes return None - elif not self.match_node_to_structure(root, structure, atoms, strict): + elif not self.match_node_to_structure(root, structure, atoms, strict, metal=metal, facet=facet, site=site): return None next_node = [] for child in root.children: - if self.match_node_to_structure(child, structure, atoms, strict): + if self.match_node_to_structure(child, structure, atoms, strict, metal=metal, facet=facet, site=site): next_node.append(child) if len(next_node) == 1: - return self.descend_tree(structure, atoms, next_node[0], strict) + return self.descend_tree(structure, atoms, next_node[0], strict, metal=metal, facet=facet, site=site) elif len(next_node) == 0: if len(root.children) > 0 and root.children[-1].label.startswith('Others-'): return root.children[-1] @@ -1072,7 +1097,7 @@ def descend_tree(self, structure, atoms, root=None, strict=False): # logging.warning('For {0}, a node {1} with overlapping children {2} was encountered ' # 'in tree with top level nodes {3}. Assuming the first match is the ' # 'better one.'.format(structure, root, next, self.top)) - return self.descend_tree(structure, atoms, next_node[0], strict) + return self.descend_tree(structure, atoms, next_node[0], strict, metal=metal, facet=facet, site=site) def are_siblings(self, node, node_other): """ @@ -1154,18 +1179,22 @@ class LogicOr(LogicNode): symbol = "OR" - def match_to_structure(self, database, structure, atoms, strict=False): + def match_to_structure(self, database, structure, atoms, strict=False, metal=None, facet=None, site=None): """ Does this node in the given database match the given structure with the labeled atoms? Setting `strict` to True makes enforces matching of atomLabels in the structure to every atomLabel in the node. + + `metal` (None or str): metal of structure (default is None) + `facet` (None or str): facet of structure (default is None) + `site` (None or str): site of structure (default is None) """ for node in self.components: if isinstance(node, LogicNode): match = node.match_to_structure(database, structure, atoms, strict) else: - match = database.match_node_to_structure(node, structure, atoms, strict) + match = database.match_node_to_structure(node, structure, atoms, strict, metal=metal, facet=facet, site=site) if match: return True != self.invert return False != self.invert @@ -1205,18 +1234,22 @@ class LogicAnd(LogicNode): symbol = "AND" - def match_to_structure(self, database, structure, atoms, strict=False): + def match_to_structure(self, database, structure, atoms, strict=False, metal=None, facet=None, site=None): """ Does this node in the given database match the given structure with the labeled atoms? Setting `strict` to True makes enforces matching of atomLabels in the structure to every atomLabel in the node. + + `metal` (None or str): metal of structure (default is None) + `facet` (None or str): facet of structure (default is None) + `site` (None or str): site of structure (default is None) """ for node in self.components: if isinstance(node, LogicNode): match = node.match_to_structure(database, structure, atoms, strict) else: - match = database.match_node_to_structure(node, structure, atoms, strict) + match = database.match_node_to_structure(node, structure, atoms, strict, metal=metal, facet=facet, site=site) if not match: return False != self.invert return True != self.invert diff --git a/rmgpy/data/baseTest.py b/rmgpy/data/baseTest.py index b1d8473d4f..b85a7e22f2 100644 --- a/rmgpy/data/baseTest.py +++ b/rmgpy/data/baseTest.py @@ -61,6 +61,17 @@ def test_match_node_to_structure(self): """) ) + entry1_facet = Entry( + item=Group().from_adjacency_list( + """ + 1 *3 C u1 {2,D} {3,S} + 2 C u0 {1,D} + 3 *5 Cd u0 {1,S} {4,D} + 4 C u0 {3,D} + """), + facet='111' + ) + entry2 = Entry( item=Group().from_adjacency_list( """ @@ -92,6 +103,17 @@ def test_match_node_to_structure(self): # entry3 contains fewer labels than entry1, therefore it can be matched self.assertTrue(self.database.match_node_to_structure(entry1, entry3.item, atoms=entry3.item.get_all_labeled_atoms())) + # entry3 does not have a facet, so it should not match entry1_facet + self.assertFalse(self.database.match_node_to_structure(entry1_facet, entry3.item, atoms=entry3.item.get_all_labeled_atoms())) + + # entry3 should match entry_1 facet if we provide a matching facet + self.assertTrue(self.database.match_node_to_structure(entry1_facet, entry3.item, atoms=entry3.item.get_all_labeled_atoms(), + facet='111')) + + # entry3 should not match entry1_facet if the facets do not match + self.assertFalse(self.database.match_node_to_structure(entry1_facet, entry3.item, atoms=entry3.item.get_all_labeled_atoms(), + facet='211')) + def test_match_node_to_node(self): """ Test that nodes can match other nodes. @@ -103,14 +125,25 @@ def test_match_node_to_node(self): """) ) + entry1_facet = Entry( + item=Group().from_adjacency_list( + """ + 1 *1 R!H u1 + """), + facet = '111' + ) + entry2 = Entry( item=Group().from_adjacency_list( """ 1 *1 Cb u1 """) ) + self.assertTrue(self.database.match_node_to_node(entry1, entry1)) + self.assertTrue(self.database.match_node_to_node(entry1_facet, entry1_facet)) self.assertFalse(self.database.match_node_to_node(entry1, entry2)) + self.assertFalse(self.database.match_node_to_node(entry1, entry1_facet)) class TestForbiddenStructures(unittest.TestCase): diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 669b6f5f6f..a4b0e12a75 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -3021,14 +3021,14 @@ def get_training_depository(self): else: raise DatabaseError('Could not find training depository in family {0}.'.format(self.label)) - def add_entry(self, parent, grp, name): + def add_entry(self, parent, grp, name, facet=None): """ Adds a group entry with parent parent group structure grp and group name name """ ind = len(self.groups.entries) - 1 - entry = Entry(index=ind, label=name, item=grp, parent=parent) + entry = Entry(index=ind, label=name, item=grp, parent=parent, facet=facet) self.groups.entries[name] = entry self.rules.entries[name] = [] if entry.parent: @@ -4401,9 +4401,9 @@ def get_reaction_matches(self, rxns=None, thermo_database=None, remove_degenerac else: mol = mol.merge(r.molecule[0]) try: - flag = not self.is_entry_match(mol, root, resonance=True) + flag = not self.is_entry_match(mol, root, resonance=True, metal=rxn.metal, facet=rxn.facet, site=rxn.site) except: - flag = not self.is_entry_match(mol, root, resonance=False) + flag = not self.is_entry_match(mol, root, resonance=False, metal=rxn.metal, facet=rxn.facet, site=rxn.site) if flag: logging.error(root.item.to_adjacency_list()) @@ -4420,7 +4420,7 @@ def get_reaction_matches(self, rxns=None, thermo_database=None, remove_degenerac while entry.children != []: for child in entry.children: - if self.is_entry_match(mol, child, resonance=False): + if self.is_entry_match(mol, child, resonance=False, metal=rxn.metal, facet=rxn.facet, site=rxn.site): entry = child rxn_lists[child.label].append(rxn) break @@ -4438,10 +4438,21 @@ def get_reaction_matches(self, rxns=None, thermo_database=None, remove_degenerac return rxn_lists - def is_entry_match(self, mol, entry, resonance=True): + def is_entry_match(self, mol, entry, resonance=True, metal=None, facet=None, site=None): """ determines if the labeled molecule object of reactants matches the entry entry """ + + if entry.metal: + if metal != entry.metal: + return False + if entry.facet: + if facet != entry.facet: + return False + if entry.site: + if site != entry.site: + return False + if isinstance(entry.item, Group): if resonance: structs = mol.generate_resonance_structures() @@ -4461,7 +4472,7 @@ def rxns_match_node(self, node, rxns): else: mol = mol.merge(r.molecule[0]) - if not self.is_entry_match(mol, node, resonance=False): + if not self.is_entry_match(mol, node, resonance=False, metal=rxn.metal, facet=rxn.facet, site=rxn.site): return False return True diff --git a/rmgpy/data/kinetics/groups.py b/rmgpy/data/kinetics/groups.py index 55c47bd263..405bc4dbd0 100644 --- a/rmgpy/data/kinetics/groups.py +++ b/rmgpy/data/kinetics/groups.py @@ -76,7 +76,7 @@ def __repr__(self): return ''.format(self.label) def load_entry(self, index, label, group, kinetics, reference=None, referenceType='', shortDesc='', longDesc='', - nodalDistance=None): + nodalDistance=None, metal=None, facet=None, site=None): """ Method for parsing entries in database files. Note that these argument names are retained for backward compatibility. @@ -103,7 +103,10 @@ def load_entry(self, index, label, group, kinetics, reference=None, referenceTyp reference_type=referenceType, short_desc=shortDesc, long_desc=longDesc.strip(), - nodal_distance=nodalDistance + nodal_distance=nodalDistance, + metal = metal, + facet = facet, + site = site ) def get_reaction_template(self, reaction): @@ -115,6 +118,7 @@ def get_reaction_template(self, reaction): # Get forward reaction template and remove any duplicates forward_template = self.top[:] + metal, facet, site = reaction.metal, reaction.facet, reaction.site temporary = [] symmetric_tree = False @@ -148,7 +152,7 @@ def get_reaction_template(self, reaction): atoms = r.get_all_labeled_atoms() - matched_node = self.descend_tree(r, atoms, root=entry, strict=True) + matched_node = self.descend_tree(r, atoms, root=entry, strict=True, metal=metal, facet=facet, site=site) if matched_node is not None: template.append(matched_node) @@ -175,7 +179,7 @@ def get_reaction_template(self, reaction): # Match structures atoms = reactant.get_all_labeled_atoms() # Descend the tree, making sure to match atomlabels exactly using strict = True - matched_node = self.descend_tree(reactant, atoms, root=entry, strict=True) + matched_node = self.descend_tree(reactant, atoms, root=entry, strict=True, metal=metal, facet=facet, site=site) if matched_node is not None: template.append(matched_node) # else: From 171b2ebb9959274736e752fd8b0bb27ad5e31aaf Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 22 Jul 2021 18:17:08 -0400 Subject: [PATCH 25/36] added facets to autotree gen extend node if we fail to split the node any furthur or all the reactions are the same, make child facet nodes --- rmgpy/data/kinetics/family.py | 34 +++++++++++++++++++++++++++++----- 1 file changed, 29 insertions(+), 5 deletions(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index a4b0e12a75..12891f5bbf 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -3269,9 +3269,12 @@ def extend_node(self, parent, template_rxn_map, obj=None, T=1000.0, iter_max=np. if exts == [] and not gave_up_split: # should only occur when all reactions at this node are identical rs = template_rxn_map[parent.label] + split_violation_flag = False for q, rxn in enumerate(rs): + if split_violation_flag: + break for j in range(q): - if not same_species_lists(rxn.reactants, rs[j].reactants, generate_initial_map=True, save_order=self.save_order): + if (not split_violation_flag) and (not same_species_lists(rxn.reactants, rs[j].reactants, generate_initial_map=True, save_order=self.save_order)): for p, atm in enumerate(parent.item.atoms): if atm.reg_dim_atm[0] != atm.reg_dim_atm[1]: logging.error('atom violation') @@ -3314,11 +3317,32 @@ def extend_node(self, parent, template_rxn_map, obj=None, T=1000.0, iter_max=np. for prod in rxn.products: logging.error(prod.label) logging.error(prod.to_adjacency_list()) - logging.error("Clearing Regularization Dimensions and Reattempting") # this usually happens when node expansion breaks some symmetry - parent.item.clear_reg_dims() # this almost always solves the problem - return True + if 'split_violations' in parent.item.props: + parent.item.props['split_violations'] += 1 + else: + parent.item.props['split_violations'] = 1 + split_violation_flag = True + if parent.item.props['split_violations'] >= 5: + # give up + # break the loop, and add facets to extend the tree + logging.error(f"Exceeded max split violations for {parent.label}") + break + else: + logging.error("Clearing Regularization Dimensions and Reattempting") # this usually happens when node expansion breaks some symmetry + parent.item.clear_reg_dims() # this almost always solves the problem + return True + # we either have all matching reactions or we broke the loop and exceeded max split violations for parent node + # lets try to extend the tree by splitting based on facet + if not parent.facet: + facets = list(set(r.facet for r in rs if r.facet)) + if len(facets) > 1: # we have more than 1 facet to split + for facet in facets: + # if (not parent.facet) or (parent.facet == facet): + # # the parent does not have a facet or the facets match, so lets create a child node with the facet + group = parent.item + extname = parent.label + '_facet{}'.format(facet) + self.add_entry(parent, group, extname, facet=facet) return False - if gave_up_split: return False From cdf047cc0ecc53c180f4cf89c026843a43441f14 Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 22 Jul 2021 18:18:46 -0400 Subject: [PATCH 26/36] do not regularize rings if the molecule contains surface site this fixed a bug where surface group parent specified rings, but children did not --- rmgpy/data/kinetics/family.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 12891f5bbf..d4755a1ec3 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -4024,7 +4024,7 @@ def simple_regularization(self, node, template_rxn_map, test=True): if not self.rxns_match_node(node, rxns): atm1.radical_electrons = oldvals - if (not skip and atm1.reg_dim_r[1] != [] and + if (not grp.contains_surface_site() and not skip and atm1.reg_dim_r[1] != [] and ('inRing' not in atm1.props.keys() or atm1.reg_dim_r[1][0] != atm1.props['inRing'])): if 'inRing' not in atm1.props.keys(): if (all(['inRing' in child.item.atoms[i].props.keys() for child in node.children]) and From 527a9d5b1a91a2455fb881de50f84c76b4f78110 Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 22 Jul 2021 18:19:30 -0400 Subject: [PATCH 27/36] added notebook to estimate A factors for surface reactions --- ipython/estimate_surface_A_factors.ipynb | 154 ++++++++++++++++++++ ipython/surface_A_factor.py | 172 +++++++++++++++++++++++ 2 files changed, 326 insertions(+) create mode 100644 ipython/estimate_surface_A_factors.ipynb create mode 100644 ipython/surface_A_factor.py diff --git a/ipython/estimate_surface_A_factors.ipynb b/ipython/estimate_surface_A_factors.ipynb new file mode 100644 index 0000000000..d134d5a13d --- /dev/null +++ b/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 +} diff --git a/ipython/surface_A_factor.py b/ipython/surface_A_factor.py new file mode 100644 index 0000000000..9f82c15315 --- /dev/null +++ b/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 From dd1eb0019b7d4a9749f8974bd238bc9be4d04e85 Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 22 Jul 2021 18:22:12 -0400 Subject: [PATCH 28/36] use cathub db branch for tests --- .github/workflows/CI.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index de5cdf8f65..a88d4ef07b 100644 --- a/.github/workflows/CI.yml +++ b/.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} From c29074628d5a19f866bf10603cbf8eaf824814c0 Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Fri, 23 Jul 2021 14:06:53 -0400 Subject: [PATCH 29/36] use `with` for mp pool --- rmgpy/data/kinetics/family.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index d4755a1ec3..4730614f1e 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -3746,9 +3746,8 @@ def make_bm_rules_from_template_rxn_map(self, template_rxn_map, nprocs=1, Tref=1 inds = inds.tolist() revinds = [inds.index(x) for x in np.arange(len(inputs))] - pool = mp.Pool(nprocs) - - kinetics_list = np.array(pool.map(_make_rule, inputs[inds])) + with mp.Pool(nprocs) as pool: + kinetics_list = np.array(pool.map(_make_rule, inputs[inds])) kinetics_list = kinetics_list[revinds] # fix order for i, kinetics in enumerate(kinetics_list): From 11c65d122232763f9acf78758860666b1851ef5c Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Fri, 23 Jul 2021 14:08:03 -0400 Subject: [PATCH 30/36] family product num is length of template products if not specified --- rmgpy/data/kinetics/family.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 4730614f1e..bbdf359100 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -736,6 +736,9 @@ def load(self, path, local_context=None, global_context=None, depository_labels= if not self.reactant_num: self.reactant_num = len(self.forward_template.reactants) + if not self.product_num: + self.product_num = len(self.forward_template.products) + self.groups.reactant_num = self.reactant_num # Generate the reverse template if necessary From 88ed85f0ceb49b07f4cf91a4d21dd01c1cb0882b Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Fri, 23 Jul 2021 14:09:34 -0400 Subject: [PATCH 31/36] raise InvalidActionError if nonexistent bond when applying recipe for `get_w0` method --- rmgpy/kinetics/arrhenius.pyx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/rmgpy/kinetics/arrhenius.pyx b/rmgpy/kinetics/arrhenius.pyx index f0229c1213..7190f13154 100644 --- a/rmgpy/kinetics/arrhenius.pyx +++ b/rmgpy/kinetics/arrhenius.pyx @@ -34,7 +34,7 @@ from scipy.optimize import curve_fit, fsolve cimport rmgpy.constants as constants import rmgpy.quantity as quantity -from rmgpy.exceptions import KineticsError +from rmgpy.exceptions import KineticsError, InvalidActionError from rmgpy.kinetics.uncertainties import rank_accuracy_map from rmgpy.kinetics import get_rate_coefficient_units_from_reaction_order from rmgpy.molecule.molecule import Bond @@ -1138,7 +1138,7 @@ def get_w0(actions, rxn): is_vdW_bond = True break if not is_vdW_bond: # no surface site, so no vdW bond - raise ('Attempted to change a nonexistent bond.') + raise InvalidActionError('Attempted to change a nonexistent bond.') else: # we found a surface site, so we will make vdw bond bd1 = Bond(a_dict[act[1]], a_dict[act[3]], order=0) else: # we have a bond From f013d37322150442a00c1f80021155b1409e86aa Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Fri, 23 Jul 2021 14:12:58 -0400 Subject: [PATCH 32/36] added `reactantNum` and `productNum` to `intra_H_migration` testing groups --- .../kinetics/families/intra_H_migration/groups.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/rmgpy/test_data/testing_database/kinetics/families/intra_H_migration/groups.py b/rmgpy/test_data/testing_database/kinetics/families/intra_H_migration/groups.py index 720f9d1443..2560cbaa8a 100644 --- a/rmgpy/test_data/testing_database/kinetics/families/intra_H_migration/groups.py +++ b/rmgpy/test_data/testing_database/kinetics/families/intra_H_migration/groups.py @@ -11,6 +11,10 @@ reversible = True +reactantNum = 1 + +productNum = 1 + autoGenerated = False recipe(actions=[ From 528c24ad0f641ef89331edc8ed11287dbf299595 Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Mon, 26 Jul 2021 12:12:24 -0400 Subject: [PATCH 33/36] saving DB metal attrs Metal attrs are loaded and saved to entries. They are loaded as Database attr, but not saved. This commit saves DB metal attrs --- rmgpy/data/base.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/rmgpy/data/base.py b/rmgpy/data/base.py index 94257f7369..6001abb5f2 100644 --- a/rmgpy/data/base.py +++ b/rmgpy/data/base.py @@ -380,6 +380,12 @@ def save(self, path, reindex=True): f.write('#!/usr/bin/env python\n') f.write('# encoding: utf-8\n\n') f.write('name = "{0}"\n'.format(self.name)) + if self.metal: + f.write('metal = "{0}"\n'.format(self.metal)) + if self.facet: + f.write('facet = "{0}"\n'.format(self.facet)) + if self.site: + f.write('site = "{0}"\n'.format(self.site)) f.write('shortDesc = "{0}"\n'.format(self.short_desc)) f.write('longDesc = """\n') f.write(self.long_desc.strip() + '\n') From 3355a3be6b2b7063d26239feea77c6887300a34c Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Mon, 26 Jul 2021 13:05:54 -0400 Subject: [PATCH 34/36] commented out `check_surface_thermo_groups_have_surface_attributes` db test The surface thermo groups should no longer have metal attrs since the adsorption groups are used for any type of metal/facet and this can mess with descending the trees if the thermo group entry has a metal, and the structure does not. Instead, the metal and facet are saved as thermo group database attrs --- testing/databaseTest.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/testing/databaseTest.py b/testing/databaseTest.py index 8be0a6e7eb..ae87ce9334 100644 --- a/testing/databaseTest.py +++ b/testing/databaseTest.py @@ -242,12 +242,12 @@ def test_thermo(self): yield test, group_name # tests for adsorption groups - if 'adsorption' in group_name.lower(): - test = lambda x: self.check_surface_thermo_groups_have_surface_attributes(group_name, group) - test_name = "Thermo surface groups {0}: Entry has metal attributes?".format(group_name) - test.description = test_name - self.compat_func_name = test_name - yield test, group_name + # if 'adsorption' in group_name.lower(): + # test = lambda x: self.check_surface_thermo_groups_have_surface_attributes(group_name, group) + # test_name = "Thermo surface groups {0}: Entry has metal attributes?".format(group_name) + # test.description = test_name + # self.compat_func_name = test_name + # yield test, group_name for library_name, library in self.database.thermo.libraries.items(): if 'surface' in library_name.lower(): From ea490709a8dd199e2da4bf1a1b475489977f51ef Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Mon, 26 Jul 2021 13:43:27 -0400 Subject: [PATCH 35/36] fixup! commented out `check_surface_thermo_groups_have_surface_attributes` db test --- testing/databaseTest.py | 48 ++++++++++++++++++++--------------------- 1 file changed, 24 insertions(+), 24 deletions(-) diff --git a/testing/databaseTest.py b/testing/databaseTest.py index ae87ce9334..82d8521ecb 100644 --- a/testing/databaseTest.py +++ b/testing/databaseTest.py @@ -1514,30 +1514,30 @@ def make_error_message(reactants, message=''): boo = True if boo: raise ValueError("Error Occurred. See log for details.") - def check_surface_thermo_groups_have_surface_attributes(self, group_name, group): - """ - Tests that each entry in the surface thermo groups has a 'metal' and 'facet' attribute, - describing which metal the data came from. - """ - failed = False - for entry in group.entries.values(): - if isinstance(entry.data, rmgpy.thermo.thermodata.ThermoData): - if 'Pt' in group_name: - if entry.metal is not 'Pt': - logging.error(f'Expected {entry} metal attribute in {group_name} group to match Pt, but was {entry.metal}') - failed = True - if '111' in group_name: - if entry.facet is not '111': - logging.error(f'Expected {entry} facet attribute in {group_name} group to match 111, but was {entry.facet}') - failed = True - if not entry.metal: - logging.error(f'Expected a metal attribute for {entry} in {group_name} group but found {entry.metal!r}') - failed = True - if not entry.facet: - logging.error(f'Expected a facet attribute for {entry} in {group_name} group but found {entry.facet!r}') - failed = True - if failed: - raise ValueError("Error occured in databaseTest. Please check log warnings for all error messages.") + # def check_surface_thermo_groups_have_surface_attributes(self, group_name, group): + # """ + # Tests that each entry in the surface thermo groups has a 'metal' and 'facet' attribute, + # describing which metal the data came from. + # """ + # failed = False + # for entry in group.entries.values(): + # if isinstance(entry.data, rmgpy.thermo.thermodata.ThermoData): + # if 'Pt' in group_name: + # if entry.metal is not 'Pt': + # logging.error(f'Expected {entry} metal attribute in {group_name} group to match Pt, but was {entry.metal}') + # failed = True + # if '111' in group_name: + # if entry.facet is not '111': + # logging.error(f'Expected {entry} facet attribute in {group_name} group to match 111, but was {entry.facet}') + # failed = True + # if not entry.metal: + # logging.error(f'Expected a metal attribute for {entry} in {group_name} group but found {entry.metal!r}') + # failed = True + # if not entry.facet: + # logging.error(f'Expected a facet attribute for {entry} in {group_name} group but found {entry.facet!r}') + # failed = True + # if failed: + # raise ValueError("Error occured in databaseTest. Please check log warnings for all error messages.") def check_surface_thermo_libraries_have_surface_attributes(self, library_name, library): """ From 5ed781e74832ead8876bdcdcc2d8e9b735efaf57 Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Mon, 26 Jul 2021 13:44:02 -0400 Subject: [PATCH 36/36] added rxn metal attrs to own reverse reactions in family `get_training_set` --- rmgpy/data/kinetics/family.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index bbdf359100..d884a9a665 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -4337,7 +4337,8 @@ def get_reactant_thermo(reactant,metal): reacts = [Species(molecule=[get_label_fixed_mol(x.molecule[0], root_labels)], thermo=x.thermo) for x in rxns[i].reactants] rrev = Reaction(reactants=products, products=reacts, - kinetics=rxns[i].generate_reverse_rate_coefficient(), rank=rxns[i].rank) + kinetics=rxns[i].generate_reverse_rate_coefficient(), rank=rxns[i].rank, + metal=rxns[i].metal, facet=rxns[i].facet, site=rxns[i].site) rrev.is_forward = False if estimate_thermo: