Merge pull request #2554 from ReactionMechanismGenerator/bidentate_Kc

fix Kc for bidentate rxns

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()

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

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

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"

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,14 +470,20 @@ 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]
         Kclist = self.rxn1s.get_equilibrium_constants(Tlist, type="Kc")
         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.