write data as RMG object instead of floats to keep the units consistent

Author: 12Chao

rmgpy/rmg/main.py

@@ -1233,8 +1233,8 @@ def execute(self, initialize=True, **kwargs):
                                 for sp in self.reaction_model.core.species:
                                     if sp.is_isomorphic(mol, strict=False):
                                         parameters['units'] = {'energy':'J', 'quantity':'mol'}
-                                        parameters['enthalpy-coefficients'] = [float(value) for value in parameters['enthalpy-coefficients']]
-                                        parameters['entropy-coefficients'] = [float(value) for value in parameters['entropy-coefficients']]
+                                        parameters['enthalpy-coefficients'] = [value.value_si for value in parameters['enthalpy-coefficients']]
+                                        parameters['entropy-coefficients'] = [value.value_si for value in parameters['entropy-coefficients']]
                                         try:
                                             content["species"][gas.n_species+surf.species_index(sp.to_chemkin())]['coverage-dependencies'][sp.to_chemkin()] = parameters
                                         except KeyError:

rmgpy/solver/surface.pyx

@@ -215,7 +215,9 @@ cdef class SurfaceReactor(ReactionSystem):
                         for species in self.species_index.keys():
                             if species.is_isomorphic(molecule, strict=False):
                                 species_index = self.species_index[species]
-                                thermo_polynomials = np.concatenate((parameters['enthalpy-coefficients'], parameters['entropy-coefficients']), axis=0)
+                                enthalpy_coeff = np.array([p.value_si for p in parameters['enthalpy-coefficients']])
+                                entropy_coeff = np.array([p.value_si for p in parameters['entropy-coefficients']])
+                                thermo_polynomials = np.concatenate((enthalpy_coeff, entropy_coeff), axis=0)
                                 self.thermo_coeff_matrix[sp_index, species_index] = [x for x in thermo_polynomials]
         # create a stoichiometry matrix for reaction enthalpy and entropy correction 
         # due to thermodynamic coverage dependence

rmgpy/thermo/nasa.pyx

@@ -35,6 +35,7 @@ from libc.math cimport log
 
 cimport rmgpy.constants as constants
 from rmgpy.util import np_list
+import rmgpy.quantity as quantity
 
 ################################################################################
 
@@ -314,8 +315,8 @@ cdef class NASA(HeatCapacityModel):
                  for species, parameters in value.items():
                     # just the polynomial model for now
                      processed_parameters = {'model': parameters['model'],
-                                             'enthalpy-coefficients': np_list([p for p in parameters['enthalpy-coefficients']]),
-                                             'entropy-coefficients': np_list([p for p in parameters['entropy-coefficients']]),
+                                             'enthalpy-coefficients': np_list([quantity.Enthalpy(p) for p in parameters['enthalpy-coefficients']]),
+                                             'entropy-coefficients': np_list([quantity.Entropy(p) for p in parameters['entropy-coefficients']]),
                                              }
                      self._thermo_coverage_dependence[species] = processed_parameters
 

rmgpy/thermo/thermodata.pyx

@@ -129,8 +129,8 @@ cdef class ThermoData(HeatCapacityModel):
                  for species, parameters in value.items():
                     # just the polynomial model for now
                      processed_parameters = {'model': parameters['model'],
-                                             'enthalpy-coefficients': np_list([p for p in parameters['enthalpy-coefficients']]),
-                                             'entropy-coefficients': np_list([p for p in parameters['entropy-coefficients']]),
+                                             'enthalpy-coefficients': np_list([quantity.Enthalpy(p) for p in parameters['enthalpy-coefficients']]),
+                                             'entropy-coefficients': np_list([quantity.Entropy(p) for p in parameters['entropy-coefficients']]),
                                              }
                      self._thermo_coverage_dependence[species] = processed_parameters
 

rmgpy/thermo/wilhoit.pyx

@@ -151,8 +151,8 @@ cdef class Wilhoit(HeatCapacityModel):
                  for species, parameters in value.items():
                     # just the polynomial model for now
                      processed_parameters = {'model': parameters['model'],
-                                             'enthalpy-coefficients': np_list([p for p in parameters['enthalpy-coefficients']]),
-                                             'entropy-coefficients': np_list([p for p in parameters['entropy-coefficients']]),
+                                             'enthalpy-coefficients': np_list([quantity.Enthalpy(p) for p in parameters['enthalpy-coefficients']]),
+                                             'entropy-coefficients': np_list([quantity.Entropy(p) for p in parameters['entropy-coefficients']]),
                                              }
                      self._thermo_coverage_dependence[species] = processed_parameters
 

test/rmgpy/thermo/nasaTest.py

@@ -31,8 +31,7 @@
 This script contains unit tests of the :mod:`rmgpy.thermo.nasa` module.
 """
 
-import os.path
-
+import os.path, ast
 
 import numpy as np
 
@@ -69,7 +68,7 @@ def setup_class(self):
         self.Tmax = 3000.0
         self.Tint = 650.73
         self.E0 = -782292.0  # J/mol.
-        self.thermo_coverage_dependence = {'1 O u0 p2 c0 {2,D} \n 2 X u0 p0 c0 {1,D}':{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,2,3]}}
+        self.thermo_coverage_dependence = {'1 O u0 p2 c0 {2,D} \n 2 X u0 p0 c0 {1,D}':{'model':'polynomial', 'enthalpy-coefficients':[(1,'J/mol'),(2,'J/mol'),(3,'J/mol')], "entropy-coefficients":[(1,'J/(mol*K)'),(2,'J/(mol*K)'),(3,'J/(mol*K)')]}}
         self.comment = "C2H6"
         self.nasa = NASA(
             polynomials=[
@@ -143,7 +142,7 @@ def test_thermo_coverage_dependence(self):
         """
         Test that the thermo_coverage_dependence property was properly set.
         """
-        assert repr(self.nasa.thermo_coverage_dependence) == repr(self.thermo_coverage_dependence)
+        assert ast.literal_eval(repr(self.nasa.thermo_coverage_dependence)) == ast.literal_eval(repr(self.thermo_coverage_dependence))
 
     def test_is_temperature_valid(self):
         """
@@ -271,7 +270,7 @@ def test_pickle(self):
         assert self.nasa.Tmax.units == nasa.Tmax.units
         assert self.nasa.E0.value == nasa.E0.value
         assert self.nasa.E0.units == nasa.E0.units
-        assert repr(self.nasa.thermo_coverage_dependence) == repr(nasa.thermo_coverage_dependence)
+        assert ast.literal_eval(repr(self.nasa.thermo_coverage_dependence)) == ast.literal_eval(repr(nasa.thermo_coverage_dependence))
         assert self.nasa.comment == nasa.comment
 
     def test_repr(self):
@@ -305,7 +304,7 @@ def test_repr(self):
         assert self.nasa.Tmax.units == nasa.Tmax.units
         assert self.nasa.E0.value == nasa.E0.value
         assert self.nasa.E0.units == nasa.E0.units
-        assert repr(self.nasa.thermo_coverage_dependence) == repr(nasa.thermo_coverage_dependence)
+        assert ast.literal_eval(repr(self.nasa.thermo_coverage_dependence)) == ast.literal_eval(repr(nasa.thermo_coverage_dependence))
         assert self.nasa.comment == nasa.comment
 
     def test_to_cantera(self):
@@ -379,7 +378,9 @@ def test_nasa_as_dict_full(self):
         assert nasa_dict["thermo_coverage_dependence"].keys() == self.thermo_coverage_dependence.keys()
         sp_name = list(self.thermo_coverage_dependence.keys())[0]
         assert nasa_dict['thermo_coverage_dependence'][sp_name]['model'] == self.thermo_coverage_dependence[sp_name]['model']
-        assert nasa_dict['thermo_coverage_dependence'][sp_name]['enthalpy-coefficients']['object'] == self.thermo_coverage_dependence[sp_name]['enthalpy-coefficients']
+        enthalpy_list = nasa_dict['thermo_coverage_dependence'][sp_name]['enthalpy-coefficients']['object'] 
+        # return [(str(coeff.value), str(coeff.units))for coeff in enthalpy_list], self.thermo_coverage_dependence[sp_name]['enthalpy-coefficients']
+        assert [(int(coeff.value), str(coeff.units))for coeff in enthalpy_list] == self.thermo_coverage_dependence[sp_name]['enthalpy-coefficients']
         assert nasa_dict['thermo_coverage_dependence'][sp_name]['entropy-coefficients']['object'] == self.thermo_coverage_dependence[sp_name]['entropy-coefficients']
         assert nasa_dict["comment"] == self.comment
         assert tuple(nasa_dict["polynomials"]["polynomial1"]["coeffs"]["object"]) == tuple(self.coeffs_low)