diff --git a/arkane/data/B2PLYP.LOG b/arkane/data/gaussian/B2PLYP.LOG
similarity index 100%
rename from arkane/data/B2PLYP.LOG
rename to arkane/data/gaussian/B2PLYP.LOG
diff --git a/arkane/data/UCCSDT_C_ATOM.LOG b/arkane/data/gaussian/UCCSDT_C_ATOM.LOG
similarity index 100%
rename from arkane/data/UCCSDT_C_ATOM.LOG
rename to arkane/data/gaussian/UCCSDT_C_ATOM.LOG
diff --git a/arkane/data/UCCSD_C_ATOM.LOG b/arkane/data/gaussian/UCCSD_C_ATOM.LOG
similarity index 100%
rename from arkane/data/UCCSD_C_ATOM.LOG
rename to arkane/data/gaussian/UCCSD_C_ATOM.LOG
diff --git a/arkane/data/UMP2_C_ATOM.LOG b/arkane/data/gaussian/UMP2_C_ATOM.LOG
similarity index 100%
rename from arkane/data/UMP2_C_ATOM.LOG
rename to arkane/data/gaussian/UMP2_C_ATOM.LOG
diff --git a/arkane/data/error_termination.out b/arkane/data/gaussian/error_termination.out
similarity index 100%
rename from arkane/data/error_termination.out
rename to arkane/data/gaussian/error_termination.out
diff --git a/arkane/data/hr_scan_with_freq.log b/arkane/data/gaussian/hr_scan_with_freq.log
similarity index 100%
rename from arkane/data/hr_scan_with_freq.log
rename to arkane/data/gaussian/hr_scan_with_freq.log
diff --git a/arkane/data/isobutanolQOOH_scan.log b/arkane/data/gaussian/isobutanolQOOH_scan.log
similarity index 100%
rename from arkane/data/isobutanolQOOH_scan.log
rename to arkane/data/gaussian/isobutanolQOOH_scan.log
diff --git a/arkane/data/Orca_TS_test.log b/arkane/data/orca/Orca_TS_test.log
similarity index 100%
rename from arkane/data/Orca_TS_test.log
rename to arkane/data/orca/Orca_TS_test.log
diff --git a/arkane/data/Orca_dlpno_test.log b/arkane/data/orca/Orca_dlpno_test.log
similarity index 100%
rename from arkane/data/Orca_dlpno_test.log
rename to arkane/data/orca/Orca_dlpno_test.log
diff --git a/arkane/data/Orca_opt_freq_test.log b/arkane/data/orca/Orca_opt_freq_test.log
similarity index 100%
rename from arkane/data/Orca_opt_freq_test.log
rename to arkane/data/orca/Orca_opt_freq_test.log
diff --git a/arkane/ess/__init__.py b/arkane/ess/__init__.py
new file mode 100644
index 00000000000..179769cbea3
--- /dev/null
+++ b/arkane/ess/__init__.py
@@ -0,0 +1,40 @@
+#!/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.                                                   #
+#                                                                             #
+###############################################################################
+
+"""
+Initialize imports for the ess (electronic structure software) modules
+"""
+
+from arkane.ess.gaussian import GaussianLog
+from arkane.ess.log import Log
+from arkane.ess.molpro import MolproLog
+from arkane.ess.orca import OrcaLog
+from arkane.ess.qchem import QChemLog
+from arkane.ess.terachem import TeraChemLog
diff --git a/arkane/logs/gaussian.py b/arkane/ess/gaussian.py
similarity index 99%
rename from arkane/logs/gaussian.py
rename to arkane/ess/gaussian.py
index a9606ecb7e3..6b35453ee66 100644
--- a/arkane/logs/gaussian.py
+++ b/arkane/ess/gaussian.py
@@ -43,7 +43,7 @@
 
 from arkane.common import check_conformer_energy, get_element_mass
 from arkane.exceptions import LogError
-from arkane.logs.log import Log
+from arkane.ess.log import Log
 
 
 ################################################################################
diff --git a/arkane/logs/gaussianTest.py b/arkane/ess/gaussianTest.py
similarity index 90%
rename from arkane/logs/gaussianTest.py
rename to arkane/ess/gaussianTest.py
index 34fb5d43c16..89bdff88a3a 100644
--- a/arkane/logs/gaussianTest.py
+++ b/arkane/ess/gaussianTest.py
@@ -28,7 +28,7 @@
 ###############################################################################
 
 """
-This module contains unit tests of the :mod:`arkane.logs.gaussian` module.
+This module contains unit tests of the :mod:`arkane.ess.gaussian` module.
 """
 
 import os
@@ -40,11 +40,10 @@
 from external.wip import work_in_progress
 from rmgpy.statmech import IdealGasTranslation, LinearRotor, NonlinearRotor, HarmonicOscillator, HinderedRotor
 
-from arkane.logs.gaussian import GaussianLog
+from arkane.ess.gaussian import GaussianLog
 from arkane.statmech import determine_qm_software
 from arkane.exceptions import LogError
 
-
 ################################################################################
 
 
@@ -91,11 +90,12 @@ def test_gaussian_energies(self):
         """
         test parsing double hydride, MP2, CCSD, CCSD(T) form Gaussian log
         """
-        log_doublehybrid = GaussianLog(os.path.join(os.path.dirname(__file__), 'data', 'B2PLYP.LOG'))
-        log_mp2 = GaussianLog(os.path.join(os.path.dirname(__file__), 'data', 'UMP2_C_ATOM.LOG'))
-        log_ccsd = GaussianLog(os.path.join(os.path.dirname(__file__), 'data', 'UCCSD_C_ATOM.LOG'))
-        log_ccsdt = GaussianLog(os.path.join(os.path.dirname(__file__), 'data', 'UCCSDT_C_ATOM.LOG'))
-        log_qb3 = GaussianLog(os.path.join(os.path.dirname(__file__), '../examples/arkane/species/C2H5/', 'ethyl_cbsqb3.log'))
+        log_doublehybrid = GaussianLog(os.path.join(self.data_path, 'B2PLYP.LOG'))
+        log_mp2 = GaussianLog(os.path.join(self.data_path, 'UMP2_C_ATOM.LOG'))
+        log_ccsd = GaussianLog(os.path.join(self.data_path, 'UCCSD_C_ATOM.LOG'))
+        log_ccsdt = GaussianLog(os.path.join(self.data_path, 'UCCSDT_C_ATOM.LOG'))
+        log_qb3 = GaussianLog(os.path.join(os.path.dirname(os.path.dirname(__file__)),
+                                           '../examples/arkane/species/C2H5/', 'ethyl_cbsqb3.log'))
 
         self.assertAlmostEqual(log_doublehybrid.load_energy() / constants.Na / constants.E_h, -0.40217794572194e+02,
                                delta=1e-6)
@@ -188,28 +188,28 @@ def test_load_scan_angle(self):
         """
         Ensures proper scan angle found in Gaussian scan job
         """
-        log = GaussianLog(os.path.join(os.path.dirname(__file__), 'data', 'isobutanolQOOH_scan.log'))
+        log = GaussianLog(os.path.join(self.data_path, 'isobutanolQOOH_scan.log'))
         self.assertAlmostEqual(log._load_scan_angle(), 10.0)
 
     def test_load_number_scans(self):
         """
         Ensures proper scan angle found in Gaussian scan job
         """
-        log = GaussianLog(os.path.join(os.path.dirname(__file__), 'data', 'isobutanolQOOH_scan.log'))
+        log = GaussianLog(os.path.join(self.data_path, 'isobutanolQOOH_scan.log'))
         self.assertAlmostEqual(log._load_number_scans(), 36)
 
     def test_determine_qm_software(self):
         """
         Ensures that determine_qm_software returns a GaussianLog object
         """
-        log = determine_qm_software(os.path.join(os.path.dirname(__file__), 'data', 'oxygen.log'))
+        log = determine_qm_software(os.path.join(self.data_path, 'oxygen.log'))
         self.assertIsInstance(log, GaussianLog)
 
     def test_gaussian_log_error_termination(self):
         """
         Ensures that error termination gaussian log file raises an logError
         """
-        file_path = os.path.join(os.path.dirname(__file__), 'data', 'error_termination.out')
+        file_path = os.path.join(self.data_path, 'error_termination.out')
         log = GaussianLog(file_path)
         with self.assertRaises(LogError) as log_error:
             log.load_conformer()
@@ -219,7 +219,7 @@ def test_load_scan_with_freq(self):
         """
         Ensures that the length of enegies with hr scans and freq calc is correct
         """
-        log = GaussianLog(os.path.join(os.path.dirname(__file__), 'data', 'hr_scan_with_freq.log'))
+        log = GaussianLog(os.path.join(self.data_path, 'hr_scan_with_freq.log'))
         self.assertAlmostEqual(log._load_number_scans(), 36)
         self.assertAlmostEqual(log._load_scan_angle(), 10.0)
         vlist, _ = log.load_scan_energies()
diff --git a/arkane/logs/log.py b/arkane/ess/log.py
similarity index 100%
rename from arkane/logs/log.py
rename to arkane/ess/log.py
diff --git a/arkane/logs/molpro.py b/arkane/ess/molpro.py
similarity index 99%
rename from arkane/logs/molpro.py
rename to arkane/ess/molpro.py
index 45b49f03f91..40db082e3ca 100644
--- a/arkane/logs/molpro.py
+++ b/arkane/ess/molpro.py
@@ -42,7 +42,7 @@
 
 from arkane.common import get_element_mass
 from arkane.exceptions import LogError
-from arkane.logs.log import Log
+from arkane.ess.log import Log
 
 ################################################################################
 
diff --git a/arkane/logs/molproTest.py b/arkane/ess/molproTest.py
similarity index 98%
rename from arkane/logs/molproTest.py
rename to arkane/ess/molproTest.py
index 786341f9b32..ad802e50569 100644
--- a/arkane/logs/molproTest.py
+++ b/arkane/ess/molproTest.py
@@ -28,7 +28,7 @@
 ###############################################################################
 
 """
-This module contains unit tests of the :mod:`arkane.logs.molpro` module.
+This module contains unit tests of the :mod:`arkane.ess.molpro` module.
 """
 
 import unittest
@@ -39,7 +39,7 @@
 import rmgpy.constants as constants
 from rmgpy.statmech import IdealGasTranslation, NonlinearRotor, HarmonicOscillator, HinderedRotor
 
-from arkane.logs.molpro import MolproLog
+from arkane.ess.molpro import MolproLog
 
 ################################################################################
 
diff --git a/arkane/logs/orca.py b/arkane/ess/orca.py
similarity index 99%
rename from arkane/logs/orca.py
rename to arkane/ess/orca.py
index d901b6476cc..5b70df0c7ab 100644
--- a/arkane/logs/orca.py
+++ b/arkane/ess/orca.py
@@ -38,7 +38,7 @@
 import rmgpy.constants as constants
 
 from arkane.common import get_element_mass
-from arkane.logs.log import Log
+from arkane.ess.log import Log
 from arkane.exceptions import LogError
 
 ################################################################################
diff --git a/arkane/logs/orcaTest.py b/arkane/ess/orcaTest.py
similarity index 79%
rename from arkane/logs/orcaTest.py
rename to arkane/ess/orcaTest.py
index 7fb6a233bb3..e72e94c80c7 100644
--- a/arkane/logs/orcaTest.py
+++ b/arkane/ess/orcaTest.py
@@ -29,13 +29,13 @@
 ###############################################################################
 
 """
-This module contains unit tests of the :mod:`arkane.logs.orca` module.
+This module contains unit tests of the :mod:`arkane.ess.orca` module.
 """
 
 import os
 import unittest
 
-from arkane.logs.orca import OrcaLog
+from arkane.ess.orca import OrcaLog
 
 ################################################################################
 
@@ -44,15 +44,21 @@ class OrcaTest(unittest.TestCase):
     """
     Contains unit tests for the orca module, used for parsing Orca log files.
     """
+    @classmethod
+    def setUpClass(cls):
+        """
+        A method that is run before all unit tests in this class.
+        """
+        cls.data_path = os.path.join(os.path.dirname(os.path.dirname(__file__)), 'data', 'orca')
 
     def test_number_of_atoms_from_orca_log(self):
         """
         Uses a Orca log files to test that
         number of atoms can be properly read.
         """
-        log = OrcaLog(os.path.join(os.path.dirname(__file__), 'data', 'Orca_opt_freq_test.log'))
+        log = OrcaLog(os.path.join(self.data_path, 'Orca_opt_freq_test.log'))
         self.assertEqual(log.get_number_of_atoms(), 3)
-        log = OrcaLog(os.path.join(os.path.dirname(__file__), 'data', 'Orca_dlpno_test.log'))
+        log = OrcaLog(os.path.join(self.data_path, 'Orca_dlpno_test.log'))
         self.assertEqual(log.get_number_of_atoms(), 30)
 
     def test_read_coordinates_from_orca_log(self):
@@ -60,10 +66,10 @@ def test_read_coordinates_from_orca_log(self):
         Uses a Orca log files to test that
         coordinate can be properly read.
         """
-        log1 = OrcaLog(os.path.join(os.path.dirname(__file__), 'data', 'Orca_opt_freq_test.log'))
+        log1 = OrcaLog(os.path.join(self.data_path, 'Orca_opt_freq_test.log'))
         coord, number, mass = log1.load_geometry()
         self.assertEqual(len(coord), 3)
-        log2 = OrcaLog(os.path.join(os.path.dirname(__file__), 'data', 'Orca_dlpno_test.log'))
+        log2 = OrcaLog(os.path.join(self.data_path, 'Orca_dlpno_test.log'))
         coord, number, mass = log2.load_geometry()
         self.assertEqual(len(coord), 30)
 
@@ -72,11 +78,11 @@ def test_energy_from_orca_log(self):
         Uses a Orca log files to test that
         molecular energies can be properly read.
         """
-        log = OrcaLog(os.path.join(os.path.dirname(__file__), 'data', 'Orca_opt_freq_test.log'))
+        log = OrcaLog(os.path.join(self.data_path, 'Orca_opt_freq_test.log'))
         self.assertAlmostEqual(log.load_energy(), -200656222.56292167, delta=1e-3)
-        log = OrcaLog(os.path.join(os.path.dirname(__file__), 'data', 'Orca_TS_test.log'))
+        log = OrcaLog(os.path.join(self.data_path, 'Orca_TS_test.log'))
         self.assertAlmostEqual(log.load_energy(), -88913623.10592663, delta=1e-3)
-        log = OrcaLog(os.path.join(os.path.dirname(__file__), 'data', 'Orca_dlpno_test.log'))
+        log = OrcaLog(os.path.join(self.data_path, 'Orca_dlpno_test.log'))
         self.assertAlmostEqual(log.load_energy(), -1816470909.1153, delta=1e-3)
 
     def test_load_zero_point_energy_from_orca_log(self):
@@ -84,9 +90,9 @@ def test_load_zero_point_energy_from_orca_log(self):
         Uses a Orca log files to test that
         molecular zero point_energy can be properly read.
         """
-        log = OrcaLog(os.path.join(os.path.dirname(__file__), 'data', 'Orca_opt_freq_test.log'))
+        log = OrcaLog(os.path.join(self.data_path, 'Orca_opt_freq_test.log'))
         self.assertAlmostEqual(log.load_zero_point_energy(), 55502.673180815, delta=1e-3)
-        log = OrcaLog(os.path.join(os.path.dirname(__file__), 'data', 'Orca_TS_test.log'))
+        log = OrcaLog(os.path.join(self.data_path, 'Orca_TS_test.log'))
         self.assertAlmostEqual(log.load_zero_point_energy(), 93500.08860598055, delta=1e-3)
 
     def test_load_negative_frequency_from_orca_log(self):
@@ -94,7 +100,7 @@ def test_load_negative_frequency_from_orca_log(self):
         Uses a orca log file for npropyl to test that its
         negative frequency can be properly read.
         """
-        log = OrcaLog(os.path.join(os.path.dirname(__file__), 'data', 'Orca_TS_test.log'))
+        log = OrcaLog(os.path.join(self.data_path, 'Orca_TS_test.log'))
         self.assertAlmostEqual(log.load_negative_frequency(), -503.24, delta=1e-1)
 
     def test_T1_diagnostic_from_orca_log(self):
@@ -102,7 +108,7 @@ def test_T1_diagnostic_from_orca_log(self):
         Uses a Orca log file for npropyl to test that its
         T1_diagnostic of freedom can be properly read.
         """
-        log = OrcaLog(os.path.join(os.path.dirname(__file__), 'data', 'Orca_dlpno_test.log'))
+        log = OrcaLog(os.path.join(self.data_path, 'Orca_dlpno_test.log'))
         self.assertAlmostEqual(log.get_T1_diagnostic(), 0.009872238, delta=1e-3)
 
 ################################################################################
diff --git a/arkane/logs/qchem.py b/arkane/ess/qchem.py
similarity index 99%
rename from arkane/logs/qchem.py
rename to arkane/ess/qchem.py
index 7b3c3a507f6..9a99debfd9d 100644
--- a/arkane/logs/qchem.py
+++ b/arkane/ess/qchem.py
@@ -43,7 +43,7 @@
 
 from arkane.common import check_conformer_energy, get_element_mass
 from arkane.exceptions import LogError
-from arkane.logs.log import Log
+from arkane.ess.log import Log
 
 ################################################################################
 
diff --git a/arkane/logs/qchemTest.py b/arkane/ess/qchemTest.py
similarity index 98%
rename from arkane/logs/qchemTest.py
rename to arkane/ess/qchemTest.py
index 76b051d72f3..870bb7aaac6 100644
--- a/arkane/logs/qchemTest.py
+++ b/arkane/ess/qchemTest.py
@@ -28,7 +28,7 @@
 ###############################################################################
 
 """
-This module contains unit tests of the :mod:`arkane.logs.qchem` module.
+This module contains unit tests of the :mod:`arkane.ess.qchem` module.
 """
 
 import os
@@ -36,7 +36,7 @@
 
 from rmgpy.statmech import IdealGasTranslation, LinearRotor, NonlinearRotor, HarmonicOscillator, HinderedRotor
 
-from arkane.logs.qchem import QChemLog
+from arkane.ess.qchem import QChemLog
 
 ################################################################################
 
diff --git a/arkane/logs/terachem.py b/arkane/ess/terachem.py
similarity index 99%
rename from arkane/logs/terachem.py
rename to arkane/ess/terachem.py
index 73c1aacf6cb..7c4c978ac26 100644
--- a/arkane/logs/terachem.py
+++ b/arkane/ess/terachem.py
@@ -43,7 +43,7 @@
 
 from arkane.common import check_conformer_energy, get_element_mass, symbol_by_number
 from arkane.exceptions import LogError
-from arkane.logs.log import Log
+from arkane.ess.log import Log
 
 ################################################################################
 
diff --git a/arkane/logs/terachemTest.py b/arkane/ess/terachemTest.py
similarity index 97%
rename from arkane/logs/terachemTest.py
rename to arkane/ess/terachemTest.py
index a5c15928201..b19621e92ab 100644
--- a/arkane/logs/terachemTest.py
+++ b/arkane/ess/terachemTest.py
@@ -28,7 +28,7 @@
 ###############################################################################
 
 """
-This module contains unit tests of the :mod:`arkane.logs.terachem` module.
+This module contains unit tests of the :mod:`arkane.ess.terachem` module.
 """
 
 import os
@@ -39,7 +39,7 @@
 from rmgpy.statmech.conformer import Conformer
 
 from arkane.exceptions import LogError
-from arkane.logs.terachem import TeraChemLog
+from arkane.ess.terachem import TeraChemLog
 
 ################################################################################
 
@@ -399,15 +399,15 @@ def test_load_scan_energies(self):
                            2.30963595e+03, 5.02046166e+03, 7.97285489e+03, 1.06923710e+04,
                            1.29244615e+04, 1.43422341e+04, 1.43905580e+04, 1.32047110e+04,
                            1.12088126e+04, 8.31162367e+03, 5.06568695e+03, 2.54966151e+03,
-                           8.50076205e+02, 0.00000000e+00, 3.31469351e+00]
-        expected_angles = [0., 0.13659098, 0.27318197, 0.40977295, 0.54636394, 0.68295492,
-                           0.81954591, 0.95613689, 1.09272788, 1.22931886, 1.36590985, 1.50250083,
-                           1.63909182, 1.7756828, 1.91227379, 2.04886477, 2.18545576, 2.32204674,
-                           2.45863773, 2.59522871, 2.7318197, 2.86841068, 3.00500167, 3.14159265,
-                           3.27818364, 3.41477462, 3.55136561, 3.68795659, 3.82454758, 3.96113856,
-                           4.09772955, 4.23432053, 4.37091152, 4.5075025, 4.64409349, 4.78068447,
-                           4.91727546, 5.05386644, 5.19045743, 5.32704841, 5.4636394, 5.60023038,
-                           5.73682137, 5.87341235, 6.01000334, 6.14659432, 6.28318531]  # radians
+                           8.50076205e+02, 0.00000000e+00]
+        expected_angles = [0., 0.13962634, 0.27925268, 0.41887902, 0.55850536, 0.6981317,
+                           0.83775804, 0.97738438, 1.11701072, 1.25663706, 1.3962634, 1.53588974,
+                           1.67551608, 1.81514242, 1.95476876, 2.0943951, 2.23402144, 2.37364778,
+                           2.51327412, 2.65290046, 2.7925268, 2.93215314, 3.07177948, 3.21140582,
+                           3.35103216, 3.4906585, 3.63028484, 3.76991118, 3.90953752, 4.04916386,
+                           4.1887902, 4.32841654, 4.46804289, 4.60766923, 4.74729557, 4.88692191,
+                           5.02654825, 5.16617459, 5.30580093, 5.44542727, 5.58505361, 5.72467995,
+                           5.86430629, 6.00393263, 6.14355897, 6.28318531]  # radians
         np.testing.assert_almost_equal(v_list, expected_v_list, 4)
         np.testing.assert_almost_equal(angles, expected_angles)
 
diff --git a/arkane/outputTest.py b/arkane/outputTest.py
index 4ae55f9c9d6..e7be8a3ef36 100644
--- a/arkane/outputTest.py
+++ b/arkane/outputTest.py
@@ -29,7 +29,7 @@
 ###############################################################################
 
 """
-This module contains unit tests of the :mod:`arkane.logs.gaussian` module.
+This module contains unit tests of the :mod:`arkane.ess.gaussian` module.
 """
 
 import os
@@ -40,7 +40,7 @@
 
 import rmgpy
 
-from arkane.logs.gaussian import GaussianLog
+from arkane.ess.gaussian import GaussianLog
 from arkane.main import Arkane
 from arkane.output import prettify, get_str_xyz
 from rmgpy.species import Species
diff --git a/arkane/statmech.py b/arkane/statmech.py
index 1da152b646a..7821138dd22 100644
--- a/arkane/statmech.py
+++ b/arkane/statmech.py
@@ -53,12 +53,7 @@
 
 from arkane.common import ArkaneSpecies, symbol_by_number, get_principal_moments_of_inertia
 from arkane.encorr.corr import get_atom_correction, get_bac
-from arkane.logs.gaussian import GaussianLog
-from arkane.logs.log import Log
-from arkane.logs.molpro import MolproLog
-from arkane.logs.orca import OrcaLog
-from arkane.logs.qchem import QChemLog
-from arkane.logs.terachem import TeraChemLog
+from arkane.ess import GaussianLog, Log, MolproLog, OrcaLog, QChemLog, TeraChemLog
 from arkane.output import prettify
 from arkane.util import determine_qm_software
 
diff --git a/arkane/statmechTest.py b/arkane/statmechTest.py
index 5c7a4f95256..28471099e9e 100644
--- a/arkane/statmechTest.py
+++ b/arkane/statmechTest.py
@@ -40,7 +40,7 @@
 from rmgpy.exceptions import InputError
 
 from arkane import Arkane
-from arkane.logs.qchem import QChemLog
+from arkane.ess.qchem import QChemLog
 from arkane.statmech import StatMechJob, determine_rotor_symmetry, is_linear
 
 ################################################################################
diff --git a/arkane/thermoTest.py b/arkane/thermoTest.py
index e14ad308523..85dcf9a591b 100644
--- a/arkane/thermoTest.py
+++ b/arkane/thermoTest.py
@@ -36,7 +36,7 @@
 
 from rmgpy.species import Species
 
-from arkane.logs.gaussian import GaussianLog
+from arkane.ess.gaussian import GaussianLog
 from arkane.thermo import ThermoJob
 
 ################################################################################
diff --git a/arkane/util.py b/arkane/util.py
index 9de74d03f74..75d5eeb5d68 100644
--- a/arkane/util.py
+++ b/arkane/util.py
@@ -35,11 +35,7 @@
 
 from rmgpy.exceptions import InputError
 
-from arkane.logs.gaussian import GaussianLog
-from arkane.logs.molpro import MolproLog
-from arkane.logs.orca import OrcaLog
-from arkane.logs.qchem import QChemLog
-from arkane.logs.terachem import TeraChemLog
+from arkane.ess import GaussianLog, MolproLog, OrcaLog, QChemLog, TeraChemLog
 
 ################################################################################
 
diff --git a/arkane/utilTest.py b/arkane/utilTest.py
index bf66a3c4ba5..f8a55b8e98f 100644
--- a/arkane/utilTest.py
+++ b/arkane/utilTest.py
@@ -36,10 +36,7 @@
 
 from rmgpy.exceptions import InputError
 
-from arkane.logs.gaussian import GaussianLog
-from arkane.logs.molpro import MolproLog
-from arkane.logs.qchem import QChemLog
-from arkane.logs.terachem import TeraChemLog
+from arkane.ess import GaussianLog, MolproLog, QChemLog, TeraChemLog
 from arkane.util import determine_qm_software
 
 ################################################################################
diff --git a/rmgpy/chemkin.pyx b/rmgpy/chemkin.pyx
index 36c53b71bf1..ebc0768c47d 100644
--- a/rmgpy/chemkin.pyx
+++ b/rmgpy/chemkin.pyx
@@ -80,7 +80,7 @@ def read_thermo_entry(entry, Tmin=0, Tint=0, Tmax=0):
 
     comment = lines[0][len(species):24].strip()
     formula = {}
-    for i in [24, 29, 34, 39, 74]:
+    for i in [24, 29, 34, 39, 73]:
         element, count = lines[0][i:i + 2].strip(), lines[0][i + 2:i + 5].strip()
         if element:
             try:
@@ -103,6 +103,21 @@ def read_thermo_entry(entry, Tmin=0, Tint=0, Tmax=0):
                         raise
             if count != 0:  # Some people put garbage elements in, with zero count. Ignore these. Allow negative counts though (eg. negative one electron)
                 formula[element] = count
+
+    # Parsing for extended elemental composition syntax, adapted from Cantera ck2cti.py
+    if lines[0].rstrip().endswith('&'):
+        complines = []
+        for i in range(len(lines)-1):
+            if lines[i].rstrip().endswith('&'):
+                complines.append(lines[i+1])
+            else:
+                break
+        lines = [lines[0]] + lines[i+1:]
+        elements = ' '.join(line.rstrip('&\n') for line in complines).split()
+        formula = {}
+        for i in range(0, len(elements), 2):
+            formula[elements[i].capitalize()] = int(elements[i+1])
+
     phase = lines[0][44]
     if phase.upper() != 'G':
         logging.warning("Was expecting gas phase thermo data for {0}. Skipping thermo data.".format(species))
@@ -1497,6 +1512,15 @@ def write_thermo_entry(species, element_counts=None, verbose=True):
     if element_counts is None:
         element_counts = get_element_count(species.molecule[0])
 
+    # Sort the element_counts dictionary so that it's C, H, Al, B, Cl, D, etc.
+    # if there's any C, else Al, B, Cl, D, H, if not. This is the "Hill" system
+    # done by Molecule.get_formula
+    if 'C' in element_counts:
+        sorted_elements = sorted(element_counts, key = lambda e: {'C':'0','H':'1'}.get(e, e))
+    else:
+        sorted_elements = sorted(element_counts)
+    element_counts = {e: element_counts[e] for e in sorted_elements}
+
     string = ''
     # Write thermo comments
     if verbose:
@@ -1509,35 +1533,44 @@ def write_thermo_entry(species, element_counts=None, verbose=True):
                 else:
                     string += "! {0}\n".format(line)
 
-    # Line 1
-    string += '{0:<16}        '.format(get_species_identifier(species))
-    if len(element_counts) <= 4:
-        # Use the original Chemkin syntax for the element counts
-        for key, count in element_counts.items():
-            if isinstance(key, tuple):
-                symbol, isotope = key
-                chemkin_name = get_element(symbol, isotope=isotope).chemkin_name
-            else:
-                chemkin_name = key
-            string += '{0!s:<2}{1:>3d}'.format(chemkin_name, count)
-        string += '     ' * (4 - len(element_counts))
-    else:
-        string += '     ' * 4
-    string += 'G{0:>10.3f}{1:>10.3f}{2:>8.2f}      1'.format(thermo.polynomials[0].Tmin.value_si,
-                                                             thermo.polynomials[1].Tmax.value_si,
-                                                             thermo.polynomials[0].Tmax.value_si)
-    if len(element_counts) > 4:
-        string += '&\n'
+    # Compile element count string
+    extended_syntax = len(element_counts) > 4  # If there are more than 4 elements, use extended syntax
+    elements = []
+    for key, count in element_counts.items():
+        if isinstance(key, tuple):
+            symbol, isotope = key
+            chemkin_name = get_element(symbol, isotope=isotope).chemkin_name
+        else:
+            chemkin_name = key
+        if extended_syntax:
+            # Create a list of alternating elements and counts
+            elements.extend([chemkin_name, str(count)])
+        else:
+            # Create a list of 5-column wide formatted element counts, e.g. 'C  10'
+            elements.append('{0!s:<2}{1:>3d}'.format(chemkin_name, count))
+    if extended_syntax:
         # Use the new-style Chemkin syntax for the element counts
+        # Place all elements in space delimited format on new line
         # This will only be recognized by Chemkin 4 or later
-        for key, count in element_counts.items():
-            if isinstance(key, tuple):
-                symbol, isotope = key
-                chemkin_name = get_element(symbol, isotope=isotope).chemkin_name
-            else:
-                chemkin_name = key
-            string += '{0!s:<2}{1:>3d}'.format(chemkin_name, count)
-    string += '\n'
+        elem_1 = ' ' * 20
+        elem_2 = '&\n' + ' '.join(elements)
+    else:
+        # Use the original Chemkin syntax for the element counts
+        # Place up to 4 elements in columns 24-43 of the first line
+        # (don't use the space in columns 74-78 for the 5th element
+        #  because nobody else does and Cantera can't read it)
+        elem_1 = ''.join(elements)
+        elem_2 = ''
+
+    # Line 1
+    string += '{ident:<16}        {elem_1:<20}G{Tmin:>10.3f}{Tint:>10.3f}{Tmax:>8.2f}      1{elem_2}\n'.format(
+        ident=get_species_identifier(species),
+        elem_1=elem_1,
+        Tmin=thermo.polynomials[0].Tmin.value_si,
+        Tint=thermo.polynomials[1].Tmax.value_si,
+        Tmax=thermo.polynomials[0].Tmax.value_si,
+        elem_2=elem_2,
+    )
 
     # Line 2
     string += '{0:< 15.8E}{1:< 15.8E}{2:< 15.8E}{3:< 15.8E}{4:< 15.8E}    2\n'.format(thermo.polynomials[1].c0,
diff --git a/rmgpy/chemkinTest.py b/rmgpy/chemkinTest.py
index fd4b4c5d2fe..59614fd8f22 100644
--- a/rmgpy/chemkinTest.py
+++ b/rmgpy/chemkinTest.py
@@ -33,7 +33,8 @@
 
 import rmgpy
 from rmgpy.chemkin import get_species_identifier, load_chemkin_file, load_transport_file, mark_duplicate_reactions, \
-    read_kinetics_entry, read_reaction_comments, read_thermo_entry, save_chemkin_file, save_species_dictionary, save_transport_file
+    read_kinetics_entry, read_reaction_comments, read_thermo_entry, save_chemkin_file, save_species_dictionary, \
+    save_transport_file, write_thermo_entry
 from rmgpy.chemkin import _remove_line_breaks, _process_duplicate_reactions
 from rmgpy.data.kinetics import LibraryReaction
 from rmgpy.exceptions import ChemkinError
@@ -41,7 +42,7 @@
 from rmgpy.kinetics.chebyshev import Chebyshev
 from rmgpy.reaction import Reaction
 from rmgpy.species import Species
-from rmgpy.thermo import NASA
+from rmgpy.thermo import NASA, NASAPolynomial
 from rmgpy.transport import TransportData
 
 
@@ -428,6 +429,120 @@ def test_mark_duplicate_reactions(self):
         self.assertEqual(duplicate_flags, expected_flags)
 
 
+class TestThermoReadWrite(unittest.TestCase):
+
+    def setUp(self):
+        """This method is run once before each test."""
+        coeffs_low = [4.03055, -0.00214171, 4.90611e-05, -5.99027e-08, 2.38945e-11, -11257.6, 3.5613]
+        coeffs_high = [-0.307954, 0.0245269, -1.2413e-05, 3.07724e-09, -3.01467e-13, -10693, 22.628]
+        Tmin = 300.
+        Tmax = 3000.
+        Tint = 650.73
+        E0 = -782292.  # J/mol.
+        comment = "C2H6"
+        self.nasa = NASA(
+            polynomials=[
+                NASAPolynomial(coeffs=coeffs_low, Tmin=(Tmin, "K"), Tmax=(Tint, "K")),
+                NASAPolynomial(coeffs=coeffs_high, Tmin=(Tint, "K"), Tmax=(Tmax, "K")),
+            ],
+            Tmin=(Tmin, "K"),
+            Tmax=(Tmax, "K"),
+            E0=(E0, "J/mol"),
+            comment=comment,
+        )
+
+        # Chemkin entries for testing - note that the values are all the same
+        self.entry1 = """C2H6                    C   2H   6          G   300.000  3000.000  650.73      1
+-3.07954000E-01 2.45269000E-02-1.24130000E-05 3.07724000E-09-3.01467000E-13    2
+-1.06930000E+04 2.26280000E+01 4.03055000E+00-2.14171000E-03 4.90611000E-05    3
+-5.99027000E-08 2.38945000E-11-1.12576000E+04 3.56130000E+00                   4
+"""
+
+        self.entry2 = """CH3NO2X                                     G   300.000  3000.000  650.73      1&
+C 1 H 3 N 1 O 2 X 1
+-3.07954000E-01 2.45269000E-02-1.24130000E-05 3.07724000E-09-3.01467000E-13    2
+-1.06930000E+04 2.26280000E+01 4.03055000E+00-2.14171000E-03 4.90611000E-05    3
+-5.99027000E-08 2.38945000E-11-1.12576000E+04 3.56130000E+00                   4
+"""
+
+        self.entry3 = """CH3NO2SX                                    G   300.000  3000.000  650.73      1&
+C 1 H 3 N 1 O 2 S 1 X 1
+-3.07954000E-01 2.45269000E-02-1.24130000E-05 3.07724000E-09-3.01467000E-13    2
+-1.06930000E+04 2.26280000E+01 4.03055000E+00-2.14171000E-03 4.90611000E-05    3
+-5.99027000E-08 2.38945000E-11-1.12576000E+04 3.56130000E+00                   4
+"""
+
+    def test_write_thermo_block(self):
+        """Test that we can write a normal thermo block"""
+        species = Species(smiles='CC')
+        species.thermo = self.nasa
+
+        result = write_thermo_entry(species, verbose=False)
+
+        self.assertEqual(result, self.entry1)
+
+    def test_read_thermo_block(self):
+        """Test that we can read a normal thermo block"""
+        species, thermo, formula = read_thermo_entry(self.entry1)
+
+        self.assertEqual(species, 'C2H6')
+        self.assertEqual(formula, {'H': 6, 'C': 2})
+        self.assertTrue(self.nasa.is_identical_to(thermo))
+
+    def test_write_thermo_block_5_elem(self):
+        """Test that we can write a thermo block for a species with 5 elements"""
+        species = Species().from_adjacency_list("""
+1 O u0 p3 c-1 {3,S}
+2 O u0 p2 c0 {3,D}
+3 N u0 p0 c+1 {1,S} {2,D} {4,S}
+4 C u0 p0 c0 {3,S} {5,S} {6,S} {7,S}
+5 H u0 p0 c0 {4,S}
+6 H u0 p0 c0 {4,S}
+7 H u0 p0 c0 {4,S}
+8 X u0 p0 c0
+""")
+        species.thermo = self.nasa
+
+        result = write_thermo_entry(species, verbose=False)
+
+        self.assertEqual(result, self.entry2)
+
+    def test_read_thermo_block_5_elem(self):
+        """Test that we can read a thermo block with 5 elements"""
+        species, thermo, formula = read_thermo_entry(self.entry2)
+
+        self.assertEqual(species, 'CH3NO2X')
+        self.assertEqual(formula, {'X': 1, 'C': 1, 'O': 2, 'H': 3, 'N': 1})
+        self.assertTrue(self.nasa.is_identical_to(thermo))
+
+    def test_write_thermo_block_6_elem(self):
+        """Test that we can write a thermo block for a species with 6 elements"""
+        species = Species().from_adjacency_list("""
+1 O u0 p3 c-1 {2,S}
+2 N u0 p0 c+1 {1,S} {3,D} {4,S}
+3 O u0 p2 c0 {2,D}
+4 C u0 p0 c0 {2,S} {5,S} {6,S} {7,S}
+5 S u0 p2 c0 {4,S} {8,S}
+6 H u0 p0 c0 {4,S}
+7 H u0 p0 c0 {4,S}
+8 H u0 p0 c0 {5,S}
+9 X u0 p0 c0
+""")
+        species.thermo = self.nasa
+
+        result = write_thermo_entry(species, verbose=False)
+
+        self.assertEqual(result, self.entry3)
+
+    def test_read_thermo_block_6_elem(self):
+        """Test that we can read a thermo block with 6 elements"""
+        species, thermo, formula = read_thermo_entry(self.entry3)
+
+        self.assertEqual(species, 'CH3NO2SX')
+        self.assertEqual(formula, {'X': 1, 'C': 1, 'O': 2, 'H': 3, 'N': 1, 'S': 1})
+        self.assertTrue(self.nasa.is_identical_to(thermo))
+
+
 class TestReadReactionComments(unittest.TestCase):
     @classmethod
     def setUpClass(cls):
diff --git a/rmgpy/data/solvation.py b/rmgpy/data/solvation.py
index fc6ff99d4b9..fd97bbea5f0 100644
--- a/rmgpy/data/solvation.py
+++ b/rmgpy/data/solvation.py
@@ -798,8 +798,8 @@ def remove_h_bonding(self, saturated_struct, added_to_radicals, added_to_pairs,
         for atom in saturated_struct.atoms:
             # Iterate over heavy (non-hydrogen) atoms
             if atom.is_non_hydrogen() and atom.radical_electrons > 0:
-                for electron in range(1, atom.radical_electrons):
-                    # Get solute data for radical group    
+                for electron in range(atom.radical_electrons):
+                    # Get solute data for radical group
                     try:
                         self._add_group_solute_data(solute_data, self.groups['radical'], saturated_struct, {'*': atom})
                     except KeyError:
diff --git a/rmgpy/data/solvationTest.py b/rmgpy/data/solvationTest.py
index 2f7a11ed97c..e7771701e19 100644
--- a/rmgpy/data/solvationTest.py
+++ b/rmgpy/data/solvationTest.py
@@ -193,6 +193,22 @@ def test_radical_and_lone_pair_generation(self):
         solute_data = self.database.get_solute_data_from_groups(species)
         self.assertIsNotNone(solute_data)
 
+    def test_radical_solute_group(self):
+        """Test that the existing radical group is found for the radical species when using group additivity"""
+        # First check whether the radical group is found for the radical species
+        rad_species = Species(smiles='[OH]')
+        rad_solute_data = self.database.get_solute_data_from_groups(rad_species)
+        self.assertTrue('radical' in rad_solute_data.comment)
+        # Then check that the radical and its saturated species give different solvation free energies
+        saturated_struct = rad_species.molecule[0].copy(deep=True)
+        saturated_struct.saturate_radicals()
+        sat_species = Species(molecule=[saturated_struct])
+        sat_solute_data = self.database.get_solute_data_from_groups(sat_species)
+        solvent_data = self.database.get_solvent_data('water')
+        rad_solvation_correction = self.database.get_solvation_correction(rad_solute_data, solvent_data)
+        sat_solvation_correction = self.database.get_solvation_correction(sat_solute_data, solvent_data)
+        self.assertNotAlmostEqual(rad_solvation_correction.gibbs / 1000, sat_solvation_correction.gibbs / 1000)
+
     def test_correction_generation(self):
         """Test we can estimate solvation thermochemistry."""
         self.testCases = [
diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py
index 74215f18496..29022fec249 100644
--- a/rmgpy/molecule/molecule.py
+++ b/rmgpy/molecule/molecule.py
@@ -1144,7 +1144,9 @@ def get_formula(self):
             symbol = atom.element.symbol
             element_dict[symbol] = element_dict.get(symbol, 0) + 1
 
-        # Use the Hill system to generate the formula
+        # Use the Hill system to generate the formula.
+        # If you change this algorithm consider also updating 
+        # the chemkin.write_thermo_entry method
         formula = ''
 
         # Carbon and hydrogen always come first if carbon is present