Stabilize thermochemistry workflow

ThermoScreening/__main__.py

@@ -6,6 +6,8 @@
 Example:
     $ python -m ThermoScreening
 """
-from cli import main
+from .cli import main
 
-main()
+
+if __name__ == "__main__":
+    main()

ThermoScreening/calculator/dftbplus.py

@@ -1,4 +1,6 @@
 import os
+import shutil
+import subprocess
 from ase.calculators.dftb import Dftb
 from ase.io import read
 import numpy as np
@@ -9,6 +11,13 @@
 # --------------------------------------------------------------------------- #
 
 
+DEFAULT_SLAKO_DIR = BASE_PATH + "../external/slakos/3ob-3-1/"
+
+
+def _slako_dir(slako_dir=None):
+    return slako_dir or os.getenv("DFTB_PREFIX") or DEFAULT_SLAKO_DIR
+
+
 class Geoopt(Dftb):
     """
     Custom DFTB+ calculator to optimize the system with the 'GeometryOptimisation' driver (Rational).
@@ -70,7 +79,7 @@ def __init__(
         super().__init__(
             atoms=atoms,
             label=label,
-            slako_dir=BASE_PATH + "../external/slakos/3ob-3-1/",
+            slako_dir=_slako_dir(slako_dir),
             Hamiltonian_Charge=charge,
             Driver_="GeometryOptimisation",
             Driver_Optimiser="Rational {}",
@@ -143,7 +152,7 @@ def __init__(
         label="second_derivative",
         charge=0,
         delta=1.0e-4,
-        slako_dir=BASE_PATH + "../external/slakos/3ob-3-1/",
+        slako_dir=None,
         **kwargs,
     ):
         """
@@ -172,7 +181,7 @@ def __init__(
         super().__init__(
             atoms=atoms,
             label=label,
-            slako_dir=slako_dir,
+            slako_dir=_slako_dir(slako_dir),
             Hamiltonian_Charge=charge,
             Driver_="SecondDerivatives",
             Driver_Delta=delta,
@@ -293,8 +302,11 @@ def calculate(self):
         None
         """
 
-        # run the modes executable
-        os.system("modes > modes.out")
+        if shutil.which("modes") is None:
+            raise FileNotFoundError("The modes executable was not found.")
+
+        with open("modes.out", "w") as output:
+            subprocess.run(["modes"], stdout=output, check=True)
 
         return None
 

ThermoScreening/cli/thermo.py

@@ -1,8 +1,7 @@
-import numpy as np
 from argparse import ArgumentParser
 from ThermoScreening.thermo.api import execute
 import time
-from ..__version__ import __version__
+from ThermoScreening.version import __version__
 
 
 def parse_args():

ThermoScreening/thermo/api.py

@@ -7,7 +7,7 @@
 from ThermoScreening import __package_name__
 
 from .inputFileReader import InputFileReader
-from .system import System
+from .system import System, dof
 from .thermo import Thermo
 from .atoms import Atom
 from ..calculator import Geoopt, Hessian, Modes
@@ -46,25 +46,18 @@ def read_xyz(coord_file: str):
         line = line.split()
         data_N = int(line[0])
 
-        # if len(line) > 1 and line[1:] is containing only numbers
-        if len(line) > 1 and all([x.replace(".", "", 1).isdigit() for x in line[1:]]):
-
-            cell = np.array(
-                [
-                    float(line[1]),
-                    float(line[2]),
-                    float(line[3]),
-                    float(line[4]),
-                    float(line[5]),
-                    float(line[6]),
-                ]
-            )
-            pbc = True
+        if len(line) >= 7:
+            try:
+                cell = np.array([float(value) for value in line[1:7]])
+            except ValueError:
+                cell = None
+            else:
+                pbc = True
 
         line = f.readline().strip()
 
         i = 0
-        data_atoms = np.zeros(data_N, dtype=str)
+        data_atoms = np.empty(data_N, dtype=object)
         data_xyz = np.zeros((data_N, 3), dtype=float)
         while True:
             line = f.readline().strip()
@@ -103,34 +96,22 @@ def read_gen(coord_file: str):
             pbc = True
         else:
             pbc = False
-        line = f.readline().strip()
-        atomic_species = str(line.split())
-        data_xyz = np.zeros((data_N, 3))
-        data_atoms = np.array([])
-        i = 0
+        atomic_species = f.readline().strip().split()
 
-        data_atoms = np.zeros(data_N, dtype=str)
+        data_atoms = np.empty(data_N, dtype=object)
         data_xyz = np.zeros((data_N, 3), dtype=float)
-        while True:
-            line = f.readline().strip()
-            line = line.split()
-            if i < data_N:
-                break
-            index = int(line[0])
+        for i in range(data_N):
+            line = f.readline().strip().split()
             symbol_number = int(line[1])
+            data_atoms[i] = atomic_species[symbol_number - 1]
             data_xyz[i, 0] = float(line[2])
             data_xyz[i, 1] = float(line[3])
             data_xyz[i, 2] = float(line[4])
-            data_atoms = np.append(str(data_atoms), atomic_species[symbol_number - 1])
-            i += 1
         if pbc:
-            line = f.readline()
-            line = f.readline()
-            line = line.split()
-            line2 = f.readline()
-            line2 = line2.split()
-            line3 = f.readline()
-            line3 = line3.split()
+            f.readline()
+            line = f.readline().split()
+            line2 = f.readline().split()
+            line3 = f.readline().split()
             cell_vector = np.array(
                 [
                     [float(line[0]), float(line[1]), float(line[2])],
@@ -201,12 +182,8 @@ def read_coord(coord_file: str, engine: str):
             data_N, data_atoms, data_xyz, cell, pbc = read_xyz(coord_file)
             return data_N, data_atoms, data_xyz, cell, pbc
         elif coord_file.endswith(".gen"):
-            # data_N, data_atoms, data_xyz, cell_vectors, pbc = read_gen(coord_file)
-            # return data_N,data_atoms,data_xyz,cell_vectors,pbc
-            logger.error(
-                "The gen file is not tested yet.",
-                exception=TSNotImplementedError
-            )
+            data_N, data_atoms, data_xyz, cell_vectors, pbc = read_gen(coord_file)
+            return data_N, data_atoms, data_xyz, cell_vectors, pbc
         else:
             logger.error(
                 "The input file is not supported.",
@@ -415,6 +392,13 @@ def run_thermo(
     for i, symbol in enumerate(atom_symbol):
         atom_list.append(Atom(symbol=symbol, position=xyz[i, :]))
 
+    expected_dof = dof(atom_list)
+    if len(vibrational_frequencies) < expected_dof:
+        logger.error(
+            "The number of vibrational frequencies does not match with the degree of freedom.",
+            exception=TSValueError
+        )
+
     system_info = System(
         atoms=atom_list,
         electronic_energy=energy,
@@ -424,12 +408,6 @@ def run_thermo(
         charge=charge,
     )
 
-    if system_info._check_frequency_length == False:
-        logger.error(
-            "The number of vibrational frequencies does not match with the degree of freedom.",
-            exception=TSValueError
-        )
-
     thermo_setup = Thermo(
         system=system_info, temperature=temperature, pressure=pressure, engine=engine
     )
@@ -549,4 +527,4 @@ def dftbplus_thermo(
         charge=charge,
     )
 
-    return thermo
+    return thermo

ThermoScreening/thermo/inputFileReader.py

@@ -69,13 +69,20 @@ def _read(self):
         -------
         None
         """
-        self._raw_input_file = open(self._input_file, "r").readlines()
+        with open(self._input_file, "r") as input_handle:
+            self._raw_input_file = input_handle.readlines()
         self._dictionary = {}
         for line in self._raw_input_file:
-            if line[0] != "#":
-                key, value = line.split(" = ")
-                self._dictionary[key.strip()] = value.strip()
-                print(key.strip(), " = ", value.strip())
+            line = line.strip()
+            if not line or line.startswith("#"):
+                continue
+            if "=" not in line:
+                self.logger.error(
+                    f"The line '{line}' is not a valid key-value assignment.",
+                    exception=TSValueError,
+                )
+            key, value = line.split("=", maxsplit=1)
+            self._dictionary[key.strip()] = value.strip()
 
         return None
 
@@ -139,4 +146,4 @@ def _check_known_keys(self):
                     exception=TSValueError
                 )
 
-        return None
+        return None

ThermoScreening/thermo/system.py

@@ -204,7 +204,10 @@ def rotational_symmetry_number(atoms: List[Atom]) -> int:
         coord.append(atom.position)
 
     mol = Molecule(name, coord)
-    return PointGroupAnalyzer(mol).get_rotational_symmetry_number
+    symmetry_number = PointGroupAnalyzer(mol).get_rotational_symmetry_number
+    if callable(symmetry_number):
+        return symmetry_number()
+    return symmetry_number
 
 
 def spacegroup_number(atoms: List[Atom], cell: Cell) -> int:

ThermoScreening/thermo/thermo.py

@@ -251,7 +251,7 @@ def _compute_rotational_partition_function(self):
         )
 
         self._rotational_partition_function = (
-            np.pi ** (1 / 2) / self._system.rotational_symmetry_number()
+            np.pi ** (1 / 2) / self._system.rotational_symmetry_number
         ) * (
             self._temperature ** (3 / 2)
             / (np.power(self._rotational_temperature_xyz, 1 / 2))

ThermoScreening/utils/header.py

@@ -1,14 +1,16 @@
-from ..__version__ import __version__
 import sys
 
+from ThermoScreening.version import __version__
+
+
 def print_header(file: str | None = None) -> None:
     """
     A function to print the header
 
     The header is printed to standard error stream. 
     """
 
-    header = f"""
+    header = rf"""
  ________   ______  
 |        \ /      \  
  \$$$$$$$$|  $$$$$$\ 
@@ -25,4 +27,4 @@ def print_header(file: str | None = None) -> None:
     if file is None:
         print(header, file=sys.stderr)
     else:
-        print(header, file=file)
+        print(header, file=file)

ThermoScreening/version.py

@@ -0,0 +1,15 @@
+from importlib.metadata import PackageNotFoundError, version
+
+
+def get_version() -> str:
+    try:
+        from .__version__ import __version__
+    except ModuleNotFoundError:
+        try:
+            return version("ThermoScreening")
+        except PackageNotFoundError:
+            return "0+unknown"
+    return __version__
+
+
+__version__ = get_version()

pyproject.toml

@@ -21,14 +21,12 @@ classifiers = [
 ]
 
 dependencies = [
-    "numpy < 2.0",
-    "argparse",
+    "numpy >= 1.26",
     "scipy",
     "pymatgen",
     "beartype",
     "ase",
     "rdkit",
-    "wfl",
 ]
 
 [project.optional-dependencies]

pytest.ini

@@ -1,2 +1,2 @@
 [pytest]
-testpaths = tests
+testpaths = tests

tests/calculator/test_dftbplus.py

@@ -1,6 +1,7 @@
 import pytest
 
 import os
+import subprocess
 import numpy as np
 import ase.io as ase_io
 from ThermoScreening import BASE_PATH
@@ -11,14 +12,42 @@
 # --------------------------------------------------------------------------- #
 
 
+def executable_starts(command):
+    try:
+        result = subprocess.run(
+            [command, "--help"],
+            capture_output=True,
+            text=True,
+            timeout=10,
+            check=False,
+        )
+    except (FileNotFoundError, subprocess.TimeoutExpired):
+        return False
+    output = result.stdout + result.stderr
+    return "DFTB+" in output and "Library not loaded" not in output
+
+
+dftbplus_ready = executable_starts("dftb+") and executable_starts("modes")
+
+
+def test_modes_missing_executable(monkeypatch):
+    monkeypatch.setattr("ThermoScreening.calculator.dftbplus.shutil.which", lambda command: None)
+    modes = Modes.__new__(Modes)
+
+    with pytest.raises(FileNotFoundError, match="modes executable"):
+        modes.calculate()
+
+
+@pytest.mark.skipif(
+    not dftbplus_ready,
+    reason="DFTB+ executables are not installed or cannot start.",
+)
 class TestDftbplus:
 
     # Test the DFTB+ calculator
     def test_dftb(self):
-        assert os.system(
-            "which dftb+ > /dev/null") == 0, "DFTB+ is not installed."
-        assert os.system(
-            "which modes > /dev/null") == 0, "Modes is not installed."
+        assert executable_starts("dftb+")
+        assert executable_starts("modes")
 
     # Test the Geoopt class
     @pytest.mark.parametrize("example_dir", ["calculator"])