From 5b159625f212d1c97f885d6cf37a1441378d2d66 Mon Sep 17 00:00:00 2001
From: tdprice-858 <>
Date: Mon, 9 Dec 2024 11:38:21 -0800
Subject: [PATCH] Added new method to thermo class which calls all internal
 methods to calculate adsorabe thermochemistry

---
 pynta/postprocessing.py | 58 +++++++++++++++++++++++++++++++++++++++--
 1 file changed, 56 insertions(+), 2 deletions(-)

diff --git a/pynta/postprocessing.py b/pynta/postprocessing.py
index e1f459be..100c3db9 100644
--- a/pynta/postprocessing.py
+++ b/pynta/postprocessing.py
@@ -719,7 +719,7 @@ def _get_vibrational_thermo(self):
 
             return
 
-        def fit_NASA(self):
+        def _fit_NASA(self):
             '''I have not yet check the math on all this stuff from Goldsmith group'''
             R = self.R
             heat_capacity = self.Cp
@@ -846,7 +846,7 @@ def fit_NASA(self):
             return
 
         # Print output in cti format
-        def format_RMG_output(self, index):
+        def _format_RMG_output(self, index):
             line = '\n'
             line += 'entry(\n    index = %s,\n' % (index)
             line += f'    label = "{self.name}",\n'
@@ -869,6 +869,60 @@ def format_RMG_output(self, index):
 
             return
 
+        def temp_corrections(self, index=0):
+            # call the subroutine for the translational and vibrational partition functions
+            _get_translation_thermo(self)
+            _get_vibrational_thermo(self)
+
+            # Atom corrections from Thermochemistry chapter by Ruscic and Bross
+            # Page 75 of DOI = https://doi.org/10.1016/B978-0-444-64087-1.00001-2
+            h_correction = 4.234  # kJ/mol. enthalpy_H(298) - enthalpy_H(0)
+            c_correction = 1.051  # kJ/mol. enthalpy_C(298) - enthalpy_C(0)
+            n_correction = 4.335  # kJ/mol. enthalpy_N(298) - enthalpy_N(0)
+            o_correction = 4.340  # kJ/mol. enthalpy_O(298) - enthalpy_O(0)
+
+            # Still not certain why these are needed. Need to read book
+            self.heat_of_formation_correction = 0.0
+            self.heat_of_formation_correction += self.composition['H'] * h_correction
+            self.heat_of_formation_correction += self.composition['C'] * c_correction
+            self.heat_of_formation_correction += self.composition['N'] * n_correction
+            self.heat_of_formation_correction += self.composition['O'] * o_correction
+
+            # note that the partition function is the production of the individual terms,
+            # whereas the thermodynamic properties are additive
+            self.Q = self.Q_trans * self.Q_vib
+            self.S = self.S_trans + self.S_vib
+            self.dH = self.dH_trans + self.dH_vib
+            self.Cp = self.Cp_trans + self.Cv_vib  # see comments in each section regarding Cp vs Cv
+            self.heat_of_formation_298K = self.heat_of_formation_0K + self.dH[
+                0] - self.heat_of_formation_correction
+            self.H = self.heat_of_formation_298K + self.dH - self.dH[0]
+
+            print(self.heat_of_formation_298K)
+            print(self.H[0])
+            # This writes H_298, S_298 and appropriate indices of Cp to file (preparation for computing adsorption corrections)
+            g = open("thermodata_adsorbates.py", 'a+')
+            g.write(
+                '[' + str(self.name) + ', Cpdata:, ' + str(self.Cp[np.where(temperature == 300)] * 239.0057)[
+                                                           1:-1] + ', ' + str(
+                    self.Cp[np.where(temperature == 400)] * 239.0057)[1:-1] + ', ' + str(
+                    self.Cp[np.where(temperature == 500)] * 239.0057)[1:-1] + ', ' + str(
+                    self.Cp[np.where(temperature == 600)] * 239.0057)[1:-1] + ', ' + str(
+                    self.Cp[np.where(temperature == 800)] * 239.0057)[1:-1] + ', ' + str(
+                    self.Cp[np.where(temperature == 1000)] * 239.0057)[1:-1] + ', ' + str(
+                    self.Cp[np.where(temperature == 1500)] * 239.0057)[
+                                                                                          1:-1] + ', ' + ",'cal/(mol*K)', H298, " + str(
+                    self.H[0] * 0.2390057) + ", 'kcal/mol', S298, " + str(
+                    self.S[0] * 239.0057) + ", 'cal/(mol*K)']")
+            g.write('\n')
+            g.close()
+
+            # now that we've computed the thermo properties, go ahead and fit them to a NASA polynomial
+            _fit_NASA(self)
+            # format_output(molecule)
+            _format_RMG_output(self, index)
+            return
+