diff --git a/rmgpy/reactionTest.py b/rmgpy/reactionTest.py index 3d8b229a88..cc16c3b8a6 100644 --- a/rmgpy/reactionTest.py +++ b/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.