added SurfaceChargeTransfer reaction tests

rmgpy/reactionTest.py

@@ -36,10 +36,11 @@
 import cantera as ct
 import numpy
 from external.wip import work_in_progress
+import math
 
 import rmgpy.constants as constants
 from rmgpy.kinetics import Arrhenius, ArrheniusEP, MultiArrhenius, PDepArrhenius, MultiPDepArrhenius, \
-    ThirdBody, Troe, Lindemann, Chebyshev, SurfaceArrhenius, StickingCoefficient
+    ThirdBody, Troe, Lindemann, Chebyshev, SurfaceArrhenius, StickingCoefficient, SurfaceChargeTransfer
 from rmgpy.molecule import Molecule
 from rmgpy.quantity import Quantity
 from rmgpy.reaction import Reaction
@@ -51,6 +52,10 @@
 from rmgpy.statmech.vibration import HarmonicOscillator
 from rmgpy.thermo import Wilhoit, ThermoData, NASA, NASAPolynomial
 
+################################################################################
+
+def order_of_magnitude(number):
+    return math.floor(math.log(number, 10))
 
 ################################################################################
 
@@ -221,6 +226,16 @@ def setUp(self):
                                       comment="""Approximate rate""")
         )
 
+    def test_n_electrons(self):
+        """Test n_electrons property"""
+        self.assertEquals(self.rxn1s.n_electrons,0)
+        self.assertEquals(self.rxn1s.n_electrons,0)
+
+    def test_n_protons(self):
+        """Test n_protons property"""
+        self.assertEquals(self.rxn1s.n_protons,0)
+        self.assertEquals(self.rxn1s.n_protons,0)
+
     def test_is_surface_reaction_species(self):
         """Test is_surface_reaction for reaction based on Species """
         self.assertTrue(self.rxn1s.is_surface_reaction())
@@ -269,6 +284,193 @@ def test_equilibrium_constant_surface_kc(self):
         for i in range(len(Tlist)):
             self.assertAlmostEqual(Kclist[i] / Kclist0[i], 1.0, 4)
 
+class TestChargeTransferReaction(unittest.TestCase):
+    """Test charge transfer reactions"""
+
+    def setUp(self):
+        m_electron = Molecule().from_smiles("e")
+        m_proton = Molecule().from_smiles("[H+]")
+        m_x = Molecule().from_adjacency_list("1 X u0 p0")
+        m_ch2x = Molecule().from_adjacency_list(
+            """
+            1 C u0 p0 c0 {2,S} {3,S} {4,D}
+            2 H u0 p0 c0 {1,S}
+            3 H u0 p0 c0 {1,S}
+            4 X u0 p0 c0 {1,D}
+            """
+            )
+        m_ch3x = Molecule().from_adjacency_list(
+            """
+            1 C u0 p0 c0 {2,S} {3,S} {4,S} {5,S}
+            2 H u0 p0 c0 {1,S}
+            3 H u0 p0 c0 {1,S}
+            4 H u0 p0 c0 {1,S}
+            5 X u0 p0 c0 {1,S}
+            """
+            )
+
+        s_electron = Species(
+            molecule=[m_electron],
+            thermo=ThermoData(Tdata=([300, 400, 500, 600, 800, 1000, 1500, 2000],
+                                     "K"),
+                              Cpdata=([0., 0., 0., 0., 0., 0., 0., 0.], "cal/(mol*K)"),
+                              H298=(0.0, "kcal/mol"),
+                              S298=(0.0, "cal/(mol*K)")))
+
+        s_proton = Species(
+            molecule=[m_proton],
+            thermo=ThermoData(Tdata=([300, 400, 500, 600, 800, 1000, 1500, 2000],
+                                     "K"),
+                              Cpdata=([0., 0., 0., 0., 0., 0., 0., 0.], "cal/(mol*K)"),
+                              H298=(0.0, "kcal/mol"),
+                              S298=(0.0, "cal/(mol*K)")))
+        s_x = Species(
+            molecule=[m_x],
+            thermo=ThermoData(Tdata=([300, 400, 500, 600, 800, 1000, 1500, 2000],
+                                     "K"),
+                              Cpdata=([0., 0., 0., 0., 0., 0., 0., 0.], "cal/(mol*K)"),
+                              H298=(0.0, "kcal/mol"),
+                              S298=(0.0, "cal/(mol*K)")))
+
+        s_ch2x = Species(
+            molecule=[m_ch2x],
+            thermo=ThermoData(Tdata=([300,400,500,600,800,1000,1500],'K'), 
+                Cpdata=([28.4959,36.3588,42.0219,46.2428,52.3978,56.921,64.1119],'J/(mol*K)'), 
+                H298=(0.654731,'kJ/mol'), S298=(19.8341,'J/(mol*K)'), Cp0=(0.01,'J/(mol*K)'), 
+                CpInf=(99.7737,'J/(mol*K)'), 
+                comment="""Thermo library: surfaceThermoPt111 Binding energy corrected by LSR (0.50C)"""))
+
+        s_ch3x = Species(
+            molecule=[m_ch3x],
+            thermo=ThermoData(Tdata=([300,400,500,600,800,1000,1500],'K'), 
+                Cpdata=([37.3325,44.9406,51.3613,56.8115,65.537,72.3287,83.3007],'J/(mol*K)'), 
+                H298=(-45.8036,'kJ/mol'), S298=(57.7449,'J/(mol*K)'), Cp0=(0.01,'J/(mol*K)'), 
+                CpInf=(124.717,'J/(mol*K)'), 
+                comment="""Thermo library: surfaceThermoPt111 Binding energy corrected by LSR (0.25C)"""))
+
+        # X=CH2 + H+ + e- <--> X-CH3
+        rxn_reduction = Reaction(reactants=[s_electron, s_proton, s_ch2x],
+                         products=[s_ch3x],
+                         kinetics=SurfaceChargeTransfer(
+                                    A = (2.483E21, 'cm^3/(mol*s)'),
+                                    V0 = (0, 'V'),
+                                    Ea = (10, 'kJ/mol'),
+                                    Tmin = (200, 'K'),
+                                    Tmax = (3000, 'K'),
+                                    ne = -1,
+                                    ))
+
+        rxn_oxidation = Reaction(products=[s_electron, s_proton, s_ch2x],
+                         reactants=[s_ch3x],
+                         kinetics=SurfaceChargeTransfer(
+                                    A = (2.483E21, 'cm^3/(mol*s)'),
+                                    V0 = (0, 'V'),
+                                    Ea = (10, 'kJ/mol'),
+                                    Tmin = (200, 'K'),
+                                    Tmax = (3000, 'K'),
+                                    ne = 1,
+                                    ))
+
+        self.rxn_reduction = rxn_reduction
+        self.rxn_oxidation = rxn_oxidation
+
+    def test_n_electrons(self):
+        """Test n_electrons property"""
+        self.assertEquals(self.rxn_reduction.n_electrons,-1)
+        self.assertEquals(self.rxn_oxidation.n_electrons,1)
+
+    def test_n_protons(self):
+        """Test n_protons property"""
+        self.assertEquals(self.rxn_reduction.n_protons,-1)
+        self.assertEquals(self.rxn_oxidation.n_protons,1)
+
+    def test_is_surface_reaction(self):
+        """Test is_surface_reaction() method"""
+        self.assertTrue(self.rxn_reduction.is_surface_reaction())
+        self.assertTrue(self.rxn_oxidation.is_surface_reaction())
+
+    def test_is_charge_transfer_reaction(self):
+        """Test is_charge_transfer_reaction() method"""
+        self.assertTrue(self.rxn_reduction.is_charge_transfer_reaction())
+        self.assertTrue(self.rxn_oxidation.is_charge_transfer_reaction())
+
+    def test_is_surface_charge_transfer(self):
+        """Test is_surface_charge_transfer() method"""
+        self.assertTrue(self.rxn_reduction.is_surface_charge_transfer_reaction())
+        self.assertTrue(self.rxn_oxidation.is_surface_charge_transfer_reaction())
+
+    def test_get_reversible_potential(self):
+        """Test get_reversible_potential() method"""
+        V0_reduction = self.rxn_reduction.get_reversible_potential(298)
+        V0_oxidation = self.rxn_oxidation.get_reversible_potential(298)
+
+        self.assertAlmostEqual(V0_reduction,V0_oxidation,6)
+        self.assertAlmostEqual(V0_reduction, 0.3967918, 6)
+        self.assertAlmostEqual(V0_oxidation, 0.3967918, 6)
+
+    def test_get_rate_coeff(self):
+        """Test get_rate_coefficient() method"""
+
+        # these should all be the same
+        kf_1 = self.rxn_reduction.get_rate_coefficient(298,0)
+        kf_2 = self.rxn_reduction.get_rate_coefficient(298,1)
+        kf_3 = self.rxn_reduction.kinetics.get_rate_coefficient(298,0)
+        self.assertAlmostEqual(kf_1, 43870506959779.0, 6)
+        self.assertAlmostEqual(kf_1, kf_2, 6)
+        self.assertAlmostEqual(kf_1, kf_3, 6)
+        self.assertAlmostEqual(kf_2, kf_3, 6)
+
+        # kf_2 should be greater than kf_1
+        kf_1 = self.rxn_oxidation.get_rate_coefficient(298,0)
+        kf_2 = self.rxn_oxidation.get_rate_coefficient(298,0.1)
+        kf_3 = self.rxn_oxidation.kinetics.get_rate_coefficient(298,0)
+        self.assertAlmostEqual(kf_1, 43870506959779.0, 6)
+        self.assertGreater(kf_2,kf_1)
+        self.assertAlmostEqual(kf_1, kf_3, 6)
+
+    def test_equilibrium_constant_surface_charge_transfer_kc(self):
+        """
+        Test the equilibrium constants of type Kc of a surface charge transfer reaction.
+        """
+        Tlist = numpy.arange(400.0, 1600.0, 200.0, numpy.float64)
+        Kclist0 = [4.19044840e+01, 3.57651471e-01, 3.16532179e-02, 7.30688731e-03,
+                    2.75080613e-03, 1.37441465e-03] #reduction
+        Kclist_reduction = self.rxn_reduction.get_equilibrium_constants(Tlist, type='Kc')
+        Kclist_oxidation = self.rxn_oxidation.get_equilibrium_constants(Tlist, type='Kc')
+        # Test a range of temperatures
+        for i in range(len(Tlist)):
+            self.assertAlmostEqual(Kclist_reduction[i] / Kclist0[i], 1.0, 6)
+            self.assertAlmostEqual(1 / Kclist_oxidation[i] / Kclist0[i], 1.0, 6)
+
+        V0 = self.rxn_oxidation.get_reversible_potential(298)
+        Kc_reduction_equil = self.rxn_reduction.get_equilibrium_constant(298, V0)
+        Kc_oxidation_equil = self.rxn_oxidation.get_equilibrium_constant(298, V0)
+        self.assertAlmostEqual(Kc_oxidation_equil * Kc_reduction_equil, 1.0, 4)
+        self.assertAlmostEqual(Kc_oxidation_equil, 1000., 4)
+        self.assertAlmostEqual(Kc_reduction_equil, 1/1000., 4)
+
+    def test_reverse_surface_charge_transfer_rate(self):
+        """
+        Test the reverse_surface_charge_transfer_rate() method
+        """
+        kf_reduction = self.rxn_reduction.kinetics
+        kf_oxidation = self.rxn_oxidation.kinetics
+        kr_oxidation = self.rxn_oxidation.generate_reverse_rate_coefficient()
+        kr_reduction = self.rxn_reduction.generate_reverse_rate_coefficient()
+        self.assertEquals(kr_reduction.A.units, 's^-1')
+        self.assertEquals(kr_oxidation.A.units, 'm^3/(mol*s)')
+
+        Tlist = numpy.linspace(298., 500., 30.,numpy.float64)
+        for T in Tlist:
+            for V in (-0.25,0.,0.25):
+                kf = kf_reduction.get_rate_coefficient(T,V)
+                kr = kr_reduction.get_rate_coefficient(T,V)
+                K = self.rxn_reduction.get_equilibrium_constant(T,V)
+                kf = kf_oxidation.get_rate_coefficient(T,V)
+                kr = kr_oxidation.get_rate_coefficient(T,V)
+                K = self.rxn_oxidation.get_equilibrium_constant(T,V)
+                self.assertEquals(order_of_magnitude(kf/kr),order_of_magnitude(K))
+
 class TestReaction(unittest.TestCase):
     """
     Contains unit tests of the Reaction class.