Merge branch 'main' into input_io_bugfix

Author: sevyharris

arkane/data/gaussian/starred_rotational_constant.log

@@ -227,6 +227,10 @@ def load_conformer(self, symmetry=None, spin_multiplicity=0, optical_isomers=Non
                     logging.debug('Conformer {0} is assigned a spin multiplicity of {1}'
                                   .format(label, spin_multiplicity))
 
+                # Keep track of rotational constants as we go in case Gaussian prints ****** instead of the value
+                if 'Rotational constants (GHZ):' in line and '*****' not in line:
+                    rot_const = [float(d) for d in line.split()[-3:]]
+
                 # The data we want is in the Thermochemistry section of the output
                 if '- Thermochemistry -' in line:
                     modes = []
@@ -246,9 +250,15 @@ def load_conformer(self, symmetry=None, spin_multiplicity=0, optical_isomers=Non
 
                         # Read moments of inertia for external rotational modes
                         elif 'Rotational constants (GHZ):' in line:
-                            inertia = [float(d) for d in line.split()[-3:]]
+                            inertia = [d for d in line.split()[-3:]]
+                            for i in range(len(inertia)):
+                                if '*******' in inertia[i]:
+                                    try:  # set rotational constant to the last known value
+                                        inertia[i] = rot_const[i]
+                                    except NameError:
+                                        raise LogError('Rotational constant listed as ***** and no previous value given')
                             for i in range(3):
-                                inertia[i] = constants.h / (8 * constants.pi * constants.pi * inertia[i] * 1e9) \
+                                inertia[i] = constants.h / (8 * constants.pi * constants.pi * float(inertia[i]) * 1e9) \
                                              * constants.Na * 1e23
                             rotation = NonlinearRotor(inertia=(inertia, "amu*angstrom^2"), symmetry=symmetry)
                             modes.append(rotation)

arkane/data/molpro/C5OH5_CCSD(T)_F12.out

@@ -350,7 +350,7 @@ def load_energy(self, zpe_scale_factor=1.):
             # Search for e_elect
             for line in lines:
                 if f12 and f12a:
-                    if 'CCSD(T)-F12a' in line and 'energy' in line:
+                    if ('CCSD(T)-F12a' in line or 'CCSD(T)-F12/' in line and '!' not in line) and 'energy' in line:
                         e_elect = float(line.split()[-1])
                         break
                 elif f12 and f12b:

arkane/data/orca/freq_orca_5.0.4.log

@@ -215,9 +215,11 @@ def load_conformer(self, symmetry=None, spin_multiplicity=0, optical_isomers=Non
             if spin_multiplicity == 0 and ' Multiplicity           Mult' in line:
                 spin_multiplicity = int(float(line.split()[3]))
 
-            if ' Mode    freq' in line:
+            if 'Mode' in line and "freq" in line:
                 frequencies = list()
                 for line_ in log[(i + 2):]:
+                    if "------" in line_:
+                        continue
                     if not line_.strip():
                         break
                     frequencies.extend([float(line_.split()[1])])
@@ -236,7 +238,7 @@ def load_conformer(self, symmetry=None, spin_multiplicity=0, optical_isomers=Non
         if len(inertia) and not all(i == 0.0 for i in inertia):
             if any(i == 0.0 for i in inertia):
                 inertia.remove(0.0)
-                rot.append(LinearRotor(inertia=(inertia, "amu*angstrom^2"), symmetry=symmetry))
+                rot.append(LinearRotor(inertia=(inertia[0], "amu*angstrom^2"), symmetry=symmetry))
             else:
                 rot.append(NonlinearRotor(inertia=(inertia, "amu*angstrom^2"), symmetry=symmetry))
 

arkane/ess/gaussian.py

@@ -105,6 +105,18 @@ dependencies:
 # configure packages to use OpenBLAS instead of Intel MKL
   - blas=*=openblas
 
+# Cantera 2.6 on the Cantera channel is configured to use Intel MKL. It doesn't specify a version minimum
+# but crashes if MKL is too old, so we have to enforce a minimum version here. But, if we can switch over
+# to the conda-forge version, we should be able to remove this line (and maybe the two after it)
+# because it use OpenBLAS instead. 
+  - conda-forge::mkl >=2023.1.0
+  # 2023.1.0 is just a guess at the minimum version. Older versions might work too.
+  - conda-forge::libcurl <= 8.9 
+  # There's a Julia/PyCall installation bug associated with libcurl 8.10+:
+  # https://discourse.julialang.org/t/curl-multi-assign-and-segmentation-fault-when-installing-package/120901/3
+  - conda-forge::pyopenssl >20
+  # ThermoCentralDatabaseInterface fails if pyopenssl is too old. 20 is just a guess at the version number. 
+
 # additional packages that are required, but not specified here (and why)
   # pydqed, pydas, mopac, and likely others require a fortran compiler (specifically gfortran)
   # in the environment. Normally we would add this to the environment file with

arkane/ess/molpro.py

@@ -62,7 +62,7 @@ def __init__(self,
                  depository=None,
                  family=None,
                  entry=None,
-                 electrons=None,
+                 electrons=0,
                  ):
         Reaction.__init__(self,
                           index=index,

arkane/ess/orca.py

@@ -246,3 +246,28 @@ def test_load_negative_frequency(self):
         with pytest.raises(LogError):
             log = GaussianLog(os.path.join(self.data_path, "rocbs-qb3_85_methanol.out"))
             imaginary_freq = log.load_negative_frequency()
+
+    def test_load_rotational_constants(self):
+        """
+        Load the rotational constants from a Gaussian output file.
+        """
+        # try reading rotational constants from a regular gaussian log file with no stars
+        log = GaussianLog(os.path.join(self.data_path, "ethylene.log"))
+        conformer, unscaled_freqs = log.load_conformer()
+        assert NonlinearRotor in [type(mode) for mode in conformer.modes]
+        for mode in conformer.modes:
+            if isinstance(mode, NonlinearRotor):
+                assert mode.symmetry == 4
+                # test the rotational constants are read in correctly
+                assert np.all(np.isclose(mode.inertia.value_si, np.array([5.69516245e-47, 2.77584017e-46, 3.34535660e-46]), atol=1e-48))
+
+
+        # test that it can read in the rotational constants even if the last instance is stars instead of digits
+        log = GaussianLog(os.path.join(self.data_path, "starred_rotational_constant.log"))
+        conformer, unscaled_freqs = log.load_conformer()
+        assert NonlinearRotor in [type(mode) for mode in conformer.modes]
+        for mode in conformer.modes:
+            if isinstance(mode, NonlinearRotor):
+                assert mode.symmetry == 1
+                # test the rotational constants are read in correctly
+                assert np.all(np.isclose(mode.inertia.value_si, np.array([5.64469824e-54, 1.86319657e-46, 1.86319657e-46]), atol=1e-48))

environment.yml

@@ -106,9 +106,12 @@ def test_load_rad_from_molpro_log_f12(self):
         """
         log = MolproLog(os.path.join(self.data_path, "OH_f12.out"))
         e0 = log.load_energy()
-
         assert round(abs(e0 / constants.Na / constants.E_h - -75.663696424380), 5) == 0
 
+        log = MolproLog(os.path.join(self.data_path, "C5OH5_CCSD(T)_F12.out"))
+        e0 = log.load_energy()
+        assert round(abs(e0 / constants.Na / constants.E_h - -268.317057640597), 5) == 0
+
     def test_load_hosi_from_molpro_log(self):
         """
         Uses a molpro log file for HOSI to test that its