diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index 3e5b4f38b4..fc3bd402a4 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -1155,6 +1155,19 @@ def contains_surface_site(self): return True return False + def number_of_surface_sites(self): + """ + Returns the number of surface sites in the molecule. + e.g. 2 for a bidentate adsorbate + """ + cython.declare(atom=Atom) + cython.declare(count=cython.int) + count = 0 + for atom in self.atoms: + if atom.is_surface_site(): + count += 1 + return count + def is_surface_site(self): """Returns ``True`` iff the molecule is nothing but a surface site 'X'.""" return len(self.atoms) == 1 and self.atoms[0].is_surface_site() diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py index ab7882cda7..886ed5c7c4 100644 --- a/rmgpy/reaction.py +++ b/rmgpy/reaction.py @@ -609,6 +609,11 @@ def get_equilibrium_constant(self, T, type='Kc', surface_site_density=2.5e-05): surface species, the `surface_site_density` is the assumed reference. """ cython.declare(dGrxn=cython.double, K=cython.double, C0=cython.double, P0=cython.double) + cython.declare(number_of_gas_reactants=cython.int, number_of_gas_products=cython.int) + cython.declare(number_of_surface_reactants=cython.int, number_of_surface_products=cython.int) + cython.declare(dN_surf=cython.int, dN_gas=cython.int, sites=cython.int) + cython.declare(sigma_nu=cython.double) + cython.declare(rectant=Species, product=Species, spcs=Species) # Use free energy of reaction to calculate Ka dGrxn = self.get_free_energy_of_reaction(T) K = np.exp(-dGrxn / constants.R / T) @@ -637,12 +642,34 @@ def get_equilibrium_constant(self, T, type='Kc', surface_site_density=2.5e-05): dN_gas = number_of_gas_products - number_of_gas_reactants # change in mols of gas spcs if type == 'Kc': + # Determine the multiplication factor of the binding site^(-stoichiometric coefficient) + # (only relevant for reactions involving multidentate adsorbates) + sigma_nu = 1 + # if there was a species with no molecule[0], then we would have presumed (above) + # that everything is gas phase, and this bit will skip. + if number_of_surface_products > 0: + for product in self.products: + sites = product.number_of_surface_sites() + if sites > 1: + # product has stoichiometric_coefficient > 0 + # so we need to divide by the number of surface sites + sigma_nu /= sites + if number_of_surface_reactants > 0: + for reactant in self.reactants: + sites = reactant.number_of_surface_sites() + if sites > 1: + # reactant has stoichiometric_coefficient < 0 + # so we need to multiply by the number of surface sites + sigma_nu *= sites + # Convert from Ka to Kc; C0 is the reference concentration if dN_gas: C0 = P0 / constants.R / T K *= C0 ** dN_gas if dN_surf: K *= surface_site_density ** dN_surf + if sigma_nu != 1: + K *= sigma_nu elif type == 'Kp': # Convert from Ka to Kp; P0 is the reference pressure K *= P0 ** dN_gas diff --git a/rmgpy/species.py b/rmgpy/species.py index 2e0df42cfe..de0c3a2d9b 100644 --- a/rmgpy/species.py +++ b/rmgpy/species.py @@ -488,6 +488,13 @@ def is_surface_site(self): """Return ``True`` if the species is a vacant surface site.""" return self.molecule[0].is_surface_site() + def number_of_surface_sites(self): + """ + Return the number of surface sites for a species. + eg. 2 for bidentate. + """ + return self.molecule[0].number_of_surface_sites() + def get_partition_function(self, T): """ Return the partition function for the species at the specified diff --git a/test/rmgpy/molecule/moleculeTest.py b/test/rmgpy/molecule/moleculeTest.py index dd21485ccd..33e75dec58 100644 --- a/test/rmgpy/molecule/moleculeTest.py +++ b/test/rmgpy/molecule/moleculeTest.py @@ -2056,6 +2056,22 @@ def test_surface_molecules(self): assert not (adsorbed.is_surface_site()) assert not (gas.is_surface_site()) + # Check the "number of surface sites" method + bidentate = Molecule().from_adjacency_list( + """ + 1 C u0 p0 c0 {2,D} {3,S} {4,S} + 2 C u0 p0 c0 {1,D} {5,S} {6,S} + 3 H u0 p0 c0 {1,S} + 4 X u0 p0 c0 {1,S} + 5 H u0 p0 c0 {2,S} + 6 X u0 p0 c0 {2,S} + """) + assert bidentate.contains_surface_site() + assert not (bidentate.is_surface_site()) + assert gas.number_of_surface_sites() == 0 + assert adsorbed.number_of_surface_sites() == 1 + assert bidentate.number_of_surface_sites() == 2 + def test_malformed_augmented_inchi(self): """Test that augmented inchi without InChI layer raises Exception.""" malform_aug_inchi = "foo" diff --git a/test/rmgpy/reactionTest.py b/test/rmgpy/reactionTest.py index 58e28eafff..26441c9e78 100644 --- a/test/rmgpy/reactionTest.py +++ b/test/rmgpy/reactionTest.py @@ -167,6 +167,16 @@ def setup_class(self): 4 H u0 p0 c0 {1,S} 5 X u0 p0 c0 {1,S}""" ) + m_cx = Molecule().from_adjacency_list( + """1 X u0 p0 c0 {2,Q} +2 C u0 p0 c0 {1,Q}""" + ) + m_cxcx = Molecule().from_adjacency_list( + """1 X u0 p0 c0 {3,D} +2 X u0 p0 c0 {4,D} +3 C u0 p0 c0 {1,D} {4,D} +4 C u0 p0 c0 {2,D} {3,D}""" + ) s_h2 = Species( label="H2(1)", @@ -284,6 +294,82 @@ def setup_class(self): ), ) + s_cx = Species( + label="CX", + molecule=[m_cx], + thermo=NASA( + polynomials=[ + NASAPolynomial( + coeffs=[ + -1.73430697E+00, + 1.89855471E-02, + -3.23563661E-05, + 2.59269890E-08, + -7.99102104E-12, + 6.36385922E+03, + 6.25445028E+00 + ], + Tmin=(298.0, 'K'), + Tmax=(1000.0, 'K') + ), + NASAPolynomial( + coeffs=[ + 2.82193241E+00, + -6.61177416E-04, + 1.24180431E-06, + -7.03993893E-10, + 1.32276605E-13, + 5.46467816E+03, + -1.55251271E+01 + ], + Tmin=(1000.0, 'K'), + Tmax=(2000.0, 'K') + ), + ], + Tmin=(298.0, 'K'), + Tmax=(2000.0, 'K'), + comment="""Thermo library: surfaceThermoPt111""", + ), + ) + + s_cxcx = Species( + label="CXCX", + molecule=[m_cxcx], + thermo=NASA( + polynomials=[ + NASAPolynomial( + coeffs=[ + 2.22282794E-01, + 1.96908600E-02, + -3.07626817E-05, + 2.35937440E-08, + -7.12438442E-12, + 2.42830208E+04, + -3.03635113E+00 + ], + Tmin=(298.0, 'K'), + Tmax=(1000.0, 'K') + ), + NASAPolynomial( + coeffs=[ + 5.60759796E+00, + -1.41921856E-03, + 2.63409811E-06, + -1.47830656E-09, + 2.75649745E-13, + 2.31084908E+04, + -2.93177601E+01 + ], + Tmin=(1000.0, 'K'), + Tmax=(2000.0, 'K') + ), + ], + Tmin=(298.0, 'K'), + Tmax=(2000.0, 'K'), + comment="""Thermo library: surfaceThermoPt111""", + ), + ) + rxn1s = Reaction( reactants=[s_h2, s_x, s_x], products=[s_hx, s_hx], @@ -339,6 +425,18 @@ def setup_class(self): ), ) + self.rxn3 = Reaction( + reactants=[s_cxcx], + products=[s_cx, s_cx], + kinetics=SurfaceArrhenius( + A=(1e13, "1/s"), + n=0.0, + Ea=(100.0, "kJ/mol"), + T0=(1.0, "K"), + comment="""Approximate rate""", + ), + ) + def test_is_surface_reaction_species(self): """Test is_surface_reaction for reaction based on Species""" assert self.rxn1s.is_surface_reaction() @@ -372,7 +470,8 @@ def test_get_rate_coefficient_units_from_reaction_order(self): def test_equilibrium_constant_surface_kc(self): """ - Test the Reaction.get_equilibrium_constant() method for Kc of a surface reaction. + Test the Reaction.get_equilibrium_constant() method for Kc of a surface reaction + and a reaction involving bidentate adsorbates """ Tlist = numpy.arange(400.0, 1600.0, 200.0, numpy.float64) Kclist0 = [15375.20186, 1.566753, 0.017772, 0.0013485, 0.000263180, 8.73504e-05] @@ -380,6 +479,11 @@ def test_equilibrium_constant_surface_kc(self): for i in range(len(Tlist)): assert round(abs(Kclist[i] / Kclist0[i] - 1.0), 4) == 0 + Kclist0_rxn3 = [8.898588e+07, 3.591402e+03, 2.279520e+01, 1.097370, 1.454379e-01, 3.436572e-02] + Kclist_rxn3 = self.rxn3.get_equilibrium_constants(Tlist, type="Kc") + for i in range(len(Tlist)): + assert round(abs(Kclist_rxn3[i] / Kclist0_rxn3[i] - 2.0), 4) == 0 + def test_reverse_sticking_coeff_rate(self): """ Test the Reaction.reverse_sticking_coeff_rate() method works for StickingCoefficient format.