Updated thermo

Author: alongd

molecule/thermo/model.pyx

@@ -163,13 +163,18 @@ cdef class HeatCapacityModel(RMGObject):
 
         Tdata = [300,400,500,600,800,1000,1500,2000]
         for T in Tdata:
-            if not (0.8 < self.get_heat_capacity(T) / other.get_heat_capacity(T) < 1.25):
+            # Do exact comparison in addition to relative in case both are zero (surface site)
+            if self.get_heat_capacity(T) != other.get_heat_capacity(T) and \
+                not (0.8 < self.get_heat_capacity(T) / other.get_heat_capacity(T) < 1.25):
                 return False
-            elif not (0.8 < self.get_enthalpy(T) / other.get_enthalpy(T) < 1.25):
+            elif self.get_enthalpy(T) != other.get_enthalpy(T) and \
+                not (0.8 < self.get_enthalpy(T) / other.get_enthalpy(T) < 1.25):
                 return False
-            elif not (0.8 < self.get_entropy(T) / other.get_entropy(T) < 1.25):
+            elif self.get_entropy(T) != other.get_entropy(T) and \
+                not (0.8 < self.get_entropy(T) / other.get_entropy(T) < 1.25):
                 return False
-            elif not (0.8 < self.get_free_energy(T) / other.get_free_energy(T) < 1.25):
+            elif self.get_free_energy(T) != other.get_free_energy(T) and \
+                not (0.8 < self.get_free_energy(T) / other.get_free_energy(T) < 1.25):
                 return False
 
         return True
@@ -185,13 +190,18 @@ cdef class HeatCapacityModel(RMGObject):
 
         Tdata = [300,400,500,600,800,1000,1500,2000]
         for T in Tdata:
-            if not (0.95 < self.get_heat_capacity(T) / other.get_heat_capacity(T) < 1.05):
+            # Do exact comparison in addition to relative in case both are zero (surface site)
+            if self.get_heat_capacity(T) != other.get_heat_capacity(T) and \
+                not (0.95 < self.get_heat_capacity(T) / other.get_heat_capacity(T) < 1.05):
                 return False
-            elif not (0.95 < self.get_enthalpy(T) / other.get_enthalpy(T) < 1.05):
+            elif self.get_enthalpy(T) != other.get_enthalpy(T) and \
+                not (0.95 < self.get_enthalpy(T) / other.get_enthalpy(T) < 1.05):
                 return False
-            elif not (0.95 < self.get_entropy(T) / other.get_entropy(T) < 1.05):
+            elif self.get_entropy(T) != other.get_entropy(T) and \
+                not (0.95 < self.get_entropy(T) / other.get_entropy(T) < 1.05):
                 return False
-            elif not (0.95 < self.get_free_energy(T) / other.get_free_energy(T) < 1.05):
+            elif self.get_free_energy(T) != other.get_free_energy(T) and \
+                not (0.95 < self.get_free_energy(T) / other.get_free_energy(T) < 1.05):
                 return False
 
         return True

molecule/thermo/thermoengine.py

@@ -37,6 +37,7 @@
 #from molecule.statmech import Conformer
 from molecule.thermo import Wilhoit, NASA, ThermoData
 from molecule.molecule import Molecule
+from molecule.molecule.fragment import Fragment
 
 
 def process_thermo_data(spc, thermo0, thermo_class=NASA, solvent_name=''):
@@ -69,7 +70,7 @@ def process_thermo_data(spc, thermo0, thermo_class=NASA, solvent_name=''):
         # correction is added to the entropy and enthalpy
         wilhoit.S0.value_si = (wilhoit.S0.value_si + solvation_correction.entropy)
         wilhoit.H0.value_si = (wilhoit.H0.value_si + solvation_correction.enthalpy)
-        wilhoit.comment += ' + Solvation correction with {} as solvent and solute estimated using {}'.format(solvent_name, solute_data.comment)
+        wilhoit.comment += f' + Solvation correction (H={solvation_correction.enthalpy/1e3:+.0f}kJ/mol;S={solvation_correction.entropy:+.0f}J/mol/K) with {solvent_name} as solvent and solute estimated using {solute_data.comment}'
 
     # Compute E0 by extrapolation to 0 K
     # if spc.conformer is None:
@@ -157,7 +158,7 @@ def evaluator(spc, solvent_name=''):
     """
     logging.debug("Evaluating spc %s ", spc)
 
-    if isinstance(spc.molecule[0], Molecule):
+    if not isinstance(spc.molecule[0], Fragment):
         spc.generate_resonance_structures()
         thermo = generate_thermo_data(spc, solvent_name=solvent_name)
     else:

molecule/thermo/wilhoit.pyx

@@ -273,8 +273,8 @@ cdef class Wilhoit(HeatCapacityModel):
 
             # What remains is to fit the polynomial coefficients (a0, a1, a2, a3)
             # This can be done directly - no iteration required
-            A = np.empty((Cpdata.shape[0],4), np.float64)
-            b = np.empty(Cpdata.shape[0], np.float64)
+            A = np.empty((Cpdata.shape[0],4), float)
+            b = np.empty(Cpdata.shape[0], float)
             for i in range(Cpdata.shape[0]):
                 y = Tdata[i] / (Tdata[i] + B)
                 for j in range(4):