Element count from conformer

arkane/commonTest.py

@@ -249,9 +249,9 @@ def testSpeciesStatmech(self):
         job.includeHinderedRotors = self.useHinderedRotors
         job.applyBondEnergyCorrections = self.useBondCorrections
         job.load()
-        self.assertTrue(isinstance(job.species.props['elementCounts'], dict))
-        self.assertEqual(job.species.props['elementCounts']['C'], 2)
-        self.assertEqual(job.species.props['elementCounts']['H'], 4)
+        self.assertTrue(isinstance(job.species.props['element_counts'], dict))
+        self.assertEqual(job.species.props['element_counts']['C'], 2)
+        self.assertEqual(job.species.props['element_counts']['H'], 4)
 
     def testSpeciesThermo(self):
         """Test thermo job execution for species from separate input file."""

arkane/statmech.py

@@ -399,7 +399,7 @@ def load(self, pdep=False):
 
         # Save atoms for use in writing thermo output
         if isinstance(self.species, Species):
-            self.species.props['elementCounts'] = atoms
+            self.species.props['element_counts'] = atoms
 
         conformer.coordinates = (coordinates, "angstroms")
         conformer.number = number

arkane/thermo.py

@@ -53,7 +53,7 @@
 from rmgpy.molecule.util import retrieveElementCount
 
 from arkane.output import prettify
-from arkane.common import ArkaneSpecies
+from arkane.common import ArkaneSpecies, symbol_by_number
 
 
 ################################################################################
@@ -155,56 +155,80 @@ def save(self, outputFile):
         species = self.species
         logging.info('Saving thermo for {0}...'.format(species.label))
 
-        f = open(outputFile, 'a')
-
-        f.write('# Thermodynamics for {0}:\n'.format(species.label))
-        H298 = species.getThermoData().getEnthalpy(298) / 4184.
-        S298 = species.getThermoData().getEntropy(298) / 4.184
-        f.write('#   Enthalpy of formation (298 K)   = {0:9.3f} kcal/mol\n'.format(H298))
-        f.write('#   Entropy of formation (298 K)    = {0:9.3f} cal/(mol*K)\n'.format(S298))
-        f.write('#    =========== =========== =========== =========== ===========\n')
-        f.write('#    Temperature Heat cap.   Enthalpy    Entropy     Free energy\n')
-        f.write('#    (K)         (cal/mol*K) (kcal/mol)  (cal/mol*K) (kcal/mol)\n')
-        f.write('#    =========== =========== =========== =========== ===========\n')
-        for T in [300, 400, 500, 600, 800, 1000, 1500, 2000, 2400]:
-            try:
-                Cp = species.getThermoData().getHeatCapacity(T) / 4.184
-                H = species.getThermoData().getEnthalpy(T) / 4184.
-                S = species.getThermoData().getEntropy(T) / 4.184
-                G = species.getThermoData().getFreeEnergy(T) / 4184.
-                f.write('#    {0:11g} {1:11.3f} {2:11.3f} {3:11.3f} {4:11.3f}\n'.format(T, Cp, H, S, G))
-            except ValueError:
-                logging.debug("Valid thermo for {0} is outside range for temperature {1}".format(species, T))
-        f.write('#    =========== =========== =========== =========== ===========\n')
-
-        thermo_string = 'thermo(label={0!r}, thermo={1!r})'.format(species.label, species.getThermoData())
-        f.write('{0}\n\n'.format(prettify(thermo_string)))
-
-        f.close()
-        # write chemkin file
-        f = open(os.path.join(os.path.dirname(outputFile), 'chem.inp'), 'a')
-        if isinstance(species, Species):
-            if species.molecule and isinstance(species.molecule[0], Molecule):
-                elementCounts = retrieveElementCount(species.molecule[0])
-            else:
+        with open(outputFile, 'a') as f:
+            f.write('# Thermodynamics for {0}:\n'.format(species.label))
+            H298 = species.getThermoData().getEnthalpy(298) / 4184.
+            S298 = species.getThermoData().getEntropy(298) / 4.184
+            f.write('#   Enthalpy of formation (298 K)   = {0:9.3f} kcal/mol\n'.format(H298))
+            f.write('#   Entropy of formation (298 K)    = {0:9.3f} cal/(mol*K)\n'.format(S298))
+            f.write('#    =========== =========== =========== =========== ===========\n')
+            f.write('#    Temperature Heat cap.   Enthalpy    Entropy     Free energy\n')
+            f.write('#    (K)         (cal/mol*K) (kcal/mol)  (cal/mol*K) (kcal/mol)\n')
+            f.write('#    =========== =========== =========== =========== ===========\n')
+            for T in [300, 400, 500, 600, 800, 1000, 1500, 2000, 2400]:
                 try:
-                    elementCounts = species.props['elementCounts']
-                except KeyError:
-                    elementCounts = {'C': 0, 'H': 0}
-        else:
-            elementCounts = {'C': 0, 'H': 0}
-        chemkin_thermo_string = writeThermoEntry(species, elementCounts=elementCounts, verbose=True)
-        f.write('{0}\n'.format(chemkin_thermo_string))
-        f.close()
+                    Cp = species.getThermoData().getHeatCapacity(T) / 4.184
+                    H = species.getThermoData().getEnthalpy(T) / 4184.
+                    S = species.getThermoData().getEntropy(T) / 4.184
+                    G = species.getThermoData().getFreeEnergy(T) / 4184.
+                    f.write('#    {0:11g} {1:11.3f} {2:11.3f} {3:11.3f} {4:11.3f}\n'.format(T, Cp, H, S, G))
+                except ValueError:
+                    logging.debug("Valid thermo for {0} is outside range for temperature {1}".format(species, T))
+            f.write('#    =========== =========== =========== =========== ===========\n')
+
+            thermo_string = 'thermo(label={0!r}, thermo={1!r})'.format(species.label, species.getThermoData())
+            f.write('{0}\n\n'.format(prettify(thermo_string)))
+
+        # write Chemkin file
+        with open(os.path.join(os.path.dirname(outputFile), 'chem.inp'), 'a') as f:
+            if isinstance(species, Species):
+                if species.molecule and isinstance(species.molecule[0], Molecule):
+                    element_counts = retrieveElementCount(species.molecule[0])
+                else:
+                    try:
+                        element_counts = species.props['element_counts']
+                    except KeyError:
+                        element_counts = self.element_count_from_conformer()
+            else:
+                element_counts = {'C': 0, 'H': 0}
+            chemkin_thermo_string = writeThermoEntry(species, elementCounts=element_counts, verbose=True)
+            f.write('{0}\n'.format(chemkin_thermo_string))
 
         # write species dictionary
         if isinstance(species, Species):
             if species.molecule and isinstance(species.molecule[0], Molecule):
-                with open(os.path.join(os.path.dirname(outputFile), 'species_dictionary.txt'), 'a') as f:
-                    f.write(species.molecule[0].toAdjacencyList(removeH=False, label=species.label))
-                    f.write('\n')
+                spec_dict_path = os.path.join(os.path.dirname(outputFile), 'species_dictionary.txt')
+                is_species_in_dict = False
+                if os.path.isfile(spec_dict_path):
+                    with open(spec_dict_path, 'r') as f:
+                        # check whether the species dictionary contains this species, in which case do not re-append
+                        for line in f.readlines():
+                            if species.label == line.strip():
+                                is_species_in_dict = True
+                                break
+                if not is_species_in_dict:
+                    with open(spec_dict_path, 'a') as f:
+                        f.write(species.molecule[0].toAdjacencyList(removeH=False, label=species.label))
+                        f.write('\n')
         return chemkin_thermo_string
 
+    def element_count_from_conformer(self):
+        """
+        Get an element count in a dictionary form (e.g.,  {'C': 3, 'H': 8}) from the species.conformer attribute.
+
+        Returns:
+            dict: Element count, keys are element symbols,
+                  values are number of occurrences of the element in the molecule.
+        """
+        element_counts = dict()
+        for number in self.species.conformer.number.value_si:
+            symbol = symbol_by_number[number]
+            if symbol in element_counts:
+                element_counts[symbol] += 1
+            else:
+                element_counts[symbol] = 1
+        return element_counts
+
     def plot(self, outputDirectory):
         """
         Plot the heat capacity, enthapy, entropy, and Gibbs free energy of the

arkane/thermoTest.py

@@ -0,0 +1,73 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+###############################################################################
+#                                                                             #
+# RMG - Reaction Mechanism Generator                                          #
+#                                                                             #
+# Copyright (c) 2002-2019 Prof. William H. Green (whgreen@mit.edu),           #
+# Prof. Richard H. West (r.west@neu.edu) and the RMG Team (rmg_dev@mit.edu)   #
+#                                                                             #
+# Permission is hereby granted, free of charge, to any person obtaining a     #
+# copy of this software and associated documentation files (the 'Software'),  #
+# to deal in the Software without restriction, including without limitation   #
+# the rights to use, copy, modify, merge, publish, distribute, sublicense,    #
+# and/or sell copies of the Software, and to permit persons to whom the       #
+# Software is furnished to do so, subject to the following conditions:        #
+#                                                                             #
+# The above copyright notice and this permission notice shall be included in  #
+# all copies or substantial portions of the Software.                         #
+#                                                                             #
+# THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR  #
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,    #
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE #
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER      #
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING     #
+# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER         #
+# DEALINGS IN THE SOFTWARE.                                                   #
+#                                                                             #
+###############################################################################
+
+"""
+This script contains unit tests of the :mod:`arkane.thermo` module.
+"""
+
+import unittest
+import os
+
+from rmgpy.species import Species
+
+from arkane.gaussian import GaussianLog
+from arkane.thermo import ThermoJob
+
+################################################################################
+
+
+class TestThermo(unittest.TestCase):
+    """
+    Contains unit tests of the ThermoJob class.
+    """
+    @classmethod
+    def setUp(cls):
+        """A method that is run before each unit test in this class"""
+        spc = Species().fromSMILES('CCO')
+        log = GaussianLog(os.path.join(os.path.dirname(__file__), 'data', 'ethylene.log'))
+        spc.conformer = log.loadConformer()[0]
+        coords, numbers, masses = log.loadGeometry()
+        spc.conformer.coordinates = coords, 'angstroms'
+        spc.conformer.number = numbers
+        spc.conformer.mass = masses, 'amu'
+        cls.thermo_job = ThermoJob(species=spc, thermoClass='NASA')
+
+    def test_element_count_from_conformer(self):
+        """Test Getting an element count dictionary from the species.conformer attribute"""
+        element_count = self.thermo_job.element_count_from_conformer()
+        self.assertEqual(element_count, {'H': 4, 'C': 2})
+
+
+
+################################################################################
+
+
+if __name__ == '__main__':
+    unittest.main(testRunner=unittest.TextTestRunner(verbosity=2))