diff --git a/ThermoScreening/calculator/dftbplus.py b/ThermoScreening/calculator/dftbplus.py
index e916ece..7c83b3f 100644
--- a/ThermoScreening/calculator/dftbplus.py
+++ b/ThermoScreening/calculator/dftbplus.py
@@ -1,15 +1,28 @@
 import os
 import shutil
 import subprocess
+
+import numpy as np
 from ase.calculators.dftb import Dftb
 from ase.io import read
-import numpy as np
 
 from ..utils.physicalConstants import PhysicalConstants
 
 # --------------------------------------------------------------------------- #
 
 
+try:
+    from PQAnalysis.analysis.vibrational import read_hessian_file as _pq_hessian_reader
+except ModuleNotFoundError:
+    _pq_hessian_reader = None
+
+
+def _read_hessian_matrix(filename):
+    if _pq_hessian_reader is not None:
+        return _pq_hessian_reader(filename)
+    return np.atleast_2d(np.loadtxt(filename, dtype=float))
+
+
 def _slako_dir(slako_dir=None):
     selected_dir = slako_dir or os.getenv("DFTB_PREFIX")
     if not selected_dir:
@@ -211,19 +224,13 @@ def read(self):
             Hessian matrix.
         """
 
-        with open("hessian.out") as f:
-            lines = [line.split() for line in f]
-
-        # matrix to array
-        hessian = []
-        for line in lines:
-            hessian += line
-        hessian = np.array(hessian, dtype=float)
+        self.hessian = _read_hessian_matrix("hessian.out")
+        hessian_size = self.atoms.get_global_number_of_atoms() * 3
 
-        self.hessian = hessian.reshape(
-            self.atoms.get_global_number_of_atoms() * 3,
-            self.atoms.get_global_number_of_atoms() * 3,
-        )
+        if self.hessian.shape != (hessian_size, hessian_size):
+            raise ValueError(
+                "Hessian matrix size does not match the number of atoms."
+            )
 
         return self.hessian
 
diff --git a/ThermoScreening/thermo/api.py b/ThermoScreening/thermo/api.py
index 538b2a1..c1b80d1 100644
--- a/ThermoScreening/thermo/api.py
+++ b/ThermoScreening/thermo/api.py
@@ -1,7 +1,9 @@
 import logging
+from pathlib import Path
 
 import numpy as np
-from PQAnalysis.io import read_gen_file
+from PQAnalysis.io import XYZFrameReader, read_gen_file
+from PQAnalysis.io.traj_file.exceptions import FrameReaderError
 
 from ThermoScreening.exceptions import TSNotImplementedError, TSValueError
 from ThermoScreening.utils.custom_logging import setup_logger
@@ -18,6 +20,23 @@
 logger = setup_logger(logger)
 
 
+def _pq_atom_names(system):
+    return np.array([atom.name for atom in system.atoms], dtype=object)
+
+
+def _pq_positions(system):
+    return np.asarray(system.pos, dtype=float)
+
+
+def _pq_xyz_cell(cell):
+    if cell.is_vacuum:
+        return None, False
+    return (
+        np.array([cell.x, cell.y, cell.z, cell.alpha, cell.beta, cell.gamma]),
+        True,
+    )
+
+
 def read_xyz(coord_file: str):
     """
     Reads xyz-file and returns number of atoms, chemical symbols of atoms, coordinates of atoms, and cell parameters.
@@ -33,45 +52,16 @@ def read_xyz(coord_file: str):
         number of atoms, chemical symbols of atoms, coordinates of atoms, cell parameters
     """
 
-    cell = None
-    pbc = False
-    data_N = None
-
-    with open(coord_file, "r") as f:
-
-        line = f.readline()
-
-        line = line.strip()
-
-        # split line and ignore spaces
-        line = line.split()
-        data_N = int(line[0])
+    try:
+        frame = XYZFrameReader().read(
+            Path(coord_file).read_text(encoding="utf-8"),
+            traj_format="xyz",
+        )
+    except (OSError, FrameReaderError, ValueError, IndexError) as exc:
+        raise TSValueError("Invalid XYZ coordinate file.") from exc
 
-        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.empty(data_N, dtype=object)
-        data_xyz = np.zeros((data_N, 3), dtype=float)
-        while True:
-            line = f.readline().strip()
-            line = line.split()
-            if len(line) == 0:
-                break
-            data_atoms[i] = str(line[0])
-            data_xyz[i, 0] = float(line[1])
-            data_xyz[i, 1] = float(line[2])
-            data_xyz[i, 2] = float(line[3])
-            i += 1
-
-    return [data_N, data_atoms, data_xyz, cell, pbc]
+    cell, pbc = _pq_xyz_cell(frame.cell)
+    return [frame.n_atoms, _pq_atom_names(frame), _pq_positions(frame), cell, pbc]
 
 
 def read_gen(coord_file: str):
@@ -91,8 +81,8 @@ def read_gen(coord_file: str):
     system = read_gen_file(coord_file)
 
     data_N = system.n_atoms
-    data_atoms = np.array([atom.name for atom in system.atoms], dtype=object)
-    data_xyz = np.asarray(system.pos, dtype=float)
+    data_atoms = _pq_atom_names(system)
+    data_xyz = _pq_positions(system)
 
     if system.cell.is_vacuum:
         cell_vector = None
@@ -120,16 +110,27 @@ def read_vib_file(vibrational_file: str):
     np.ndarray
         The vibrational frequencies as a numpy array.
     """
-    vibrational_frequencies = np.array([])
-    with open(vibrational_file, "r") as f:
-        while True:
-            line = f.readline().strip()
-            line = line.split()
-            if len(line) == 0:
-                break
-            vibrational_frequencies = np.append(vibrational_frequencies, float(line[1]))
+    vibrational_frequencies = []
+    with open(vibrational_file, "r", encoding="utf-8") as f:
+        for line_number, line in enumerate(f, start=1):
+            fields = line.split()
+            if not fields:
+                continue
+            if len(fields) < 2:
+                raise TSValueError(
+                    f"Invalid vibrational frequency line {line_number}."
+                )
+            try:
+                vibrational_frequencies.append(float(fields[1]))
+            except ValueError as exc:
+                raise TSValueError(
+                    f"Invalid vibrational frequency line {line_number}."
+                ) from exc
 
-    return vibrational_frequencies
+    if not vibrational_frequencies:
+        raise TSValueError("No vibrational frequencies found.")
+
+    return np.asarray(vibrational_frequencies)
 
 
 def read_coord(coord_file: str, engine: str):
diff --git a/ThermoScreening/thermo/atoms.py b/ThermoScreening/thermo/atoms.py
index f19781d..4f4d8ef 100644
--- a/ThermoScreening/thermo/atoms.py
+++ b/ThermoScreening/thermo/atoms.py
@@ -1,6 +1,11 @@
 import logging
 
 import numpy as np
+from PQAnalysis.core.atom.element import (
+    atomicMasses,
+    atomicNumbers,
+    atomicNumbersReverse,
+)
 
 from ThermoScreening.exceptions import TSValueError
 from ThermoScreening.utils.custom_logging import setup_logger
@@ -243,341 +248,7 @@ def change_atom(
             self.position = position        
 
 
-atomicMasses = {
-    "h": 1.00794,
-    "d": 2.014101778,
-    "t": 3.0160492675,
-    "he": 4.002602,
-    "li": 6.941,
-    "be": 9.012182,
-    "b": 10.811,
-    "c": 12.0107,
-    "n": 14.0067,
-    "o": 15.9994,
-    "f": 18.9984032,
-    "ne": 20.1797,
-    "na": 22.989770,
-    "mg": 24.3050,
-    "al": 26.981538,
-    "si": 28.0855,
-    "p": 30.973761,
-    "s": 32.065,
-    "cl": 35.453,
-    "ar": 39.948,
-    "k": 39.0983,
-    "ca": 40.078,
-    "sc": 44.955910,
-    "ti": 47.880,
-    "v": 50.9415,
-    "cr": 51.9961,
-    "mn": 54.938049,
-    "fe": 55.845,
-    "co": 58.933200,
-    "ni": 58.6934,
-    "cu": 63.546,
-    "zn": 65.399,
-    "ga": 69.723,
-    "ge": 72.64,
-    "as": 74.92160,
-    "se": 78.96,
-    "br": 79.904,
-    "kr": 83.798,
-    "rb": 85.4678,
-    "sr": 87.62,
-    "y": 88.90585,
-    "zr": 91.224,
-    "nb": 92.90638,
-    "mo": 95.94,
-    "tc": 98.9063,
-    "ru": 101.07,
-    "rh": 102.9055,
-    "pd": 106.42,
-    "ag": 107.8682,
-    "cd": 112.411,
-    "in": 114.818,
-    "sn": 118.71,
-    "sb": 121.76,
-    "te": 127.6,
-    "i": 126.90447,
-    "xe": 131.293,
-    "cs": 132.90546,
-    "ba": 137.327,
-    "la": 138.9055,
-    "ce": 140.116,
-    "pr": 140.90765,
-    "nd": 144.24,
-    "pm": 146.9151,
-    "sm": 150.36,
-    "eu": 151.964,
-    "gd": 157.25,
-    "tb": 158.92534,
-    "dy": 162.5,
-    "ho": 164.93032,
-    "er": 167.259,
-    "tm": 168.93421,
-    "yb": 173.04,
-    "lu": 174.967,
-    "hf": 178.49,
-    "ta": 180.9479,
-    "w": 183.84,
-    "re": 186.207,
-    "os": 190.23,
-    "ir": 192.217,
-    "pt": 195.078,
-    "au": 196.96655,
-    "hg": 200.59,
-    "tl": 204.3833,
-    "pb": 207.2,
-    "bi": 208.98038,
-    "po": 208.9824,
-    "at": 209.9871,
-    "rn": 222.0176,
-    "fr": 223.0197,
-    "ra": 226.0254,
-    "ac": 227.0278,
-    "th": 232.0381,
-    "pa": 231.03588,
-    "u": 238.0289,
-    "np": 237.0482,
-    "pu": 244.0642,
-    "am": 243.0614,
-    "cm": 247.0703,
-    "bk": 247.0703,
-    "cf": 251.0796,
-    "es": 252.0829,
-    "fm": 257.0951,
-    "md": 258.0986,
-    "no": 259.1009,
-    "lr": 260.1053,
-    "q": 999.00000,
-    "x": 999.00000,
-    "cav": 1000.00000,
-    "sup": 1000000.0,
-    "dum": 1.0,
-}
-
-atomicNumbers = {
-    "h": 1,
-    "d": 1,
-    "t": 1,
-    "he": 2,
-    "li": 3,
-    "be": 4,
-    "b": 5,
-    "c": 6,
-    "n": 7,
-    "o": 8,
-    "f": 9,
-    "ne": 10,
-    "na": 11,
-    "mg": 12,
-    "al": 13,
-    "si": 14,
-    "p": 15,
-    "s": 16,
-    "cl": 17,
-    "ar": 18,
-    "k": 19,
-    "ca": 20,
-    "sc": 21,
-    "ti": 22,
-    "v": 23,
-    "cr": 24,
-    "mn": 25,
-    "fe": 26,
-    "co": 27,
-    "ni": 28,
-    "cu": 29,
-    "zn": 30,
-    "ga": 31,
-    "ge": 32,
-    "as": 33,
-    "se": 34,
-    "br": 35,
-    "kr": 36,
-    "rb": 37,
-    "sr": 38,
-    "y": 39,
-    "zr": 40,
-    "nb": 41,
-    "mo": 42,
-    "tc": 43,
-    "ru": 44,
-    "rh": 45,
-    "pd": 46,
-    "ag": 47,
-    "cd": 48,
-    "in": 49,
-    "sn": 50,
-    "sb": 51,
-    "te": 52,
-    "i": 53,
-    "xe": 54,
-    "cs": 55,
-    "ba": 56,
-    "la": 57,
-    "ce": 58,
-    "pr": 59,
-    "nd": 60,
-    "pm": 61,
-    "sm": 62,
-    "eu": 63,
-    "gd": 64,
-    "tb": 65,
-    "dy": 66,
-    "ho": 67,
-    "er": 68,
-    "tm": 69,
-    "yb": 70,
-    "lu": 71,
-    "hf": 72,
-    "ta": 73,
-    "w": 74,
-    "re": 75,
-    "os": 76,
-    "ir": 77,
-    "pt": 78,
-    "au": 79,
-    "hg": 80,
-    "tl": 81,
-    "pb": 82,
-    "bi": 83,
-    "po": 84,
-    "at": 85,
-    "rn": 86,
-    "fr": 87,
-    "ra": 88,
-    "ac": 89,
-    "th": 90,
-    "pa": 91,
-    "u": 92,
-    "np": 93,
-    "pu": 94,
-    "am": 95,
-    "cm": 96,
-    "bk": 97,
-    "cf": 98,
-    "es": 99,
-    "fm": 100,
-    "md": 101,
-    "no": 102,
-    "lr": 103,
-    "q": 999,
-    "x": 999,
-    "cav": 1000,
-    "sup": 1000000,
-    "dum": 1,
-}
-atomic_Symbol = {
-    1: "h",
-    2: "he",
-    3: "li",
-    4: "be",
-    5: "b",
-    6: "c",
-    7: "n",
-    8: "o",
-    9: "f",
-    10: "ne",
-    11: "na",
-    12: "mg",
-    13: "al",
-    14: "si",
-    15: "p",
-    16: "s",
-    17: "cl",
-    18: "ar",
-    19: "k",
-    20: "ca",
-    21: "sc",
-    22: "ti",
-    23: "v",
-    24: "cr",
-    25: "mn",
-    26: "fe",
-    27: "co",
-    28: "ni",
-    29: "cu",
-    30: "zn",
-    31: "ga",
-    32: "ge",
-    33: "as",
-    34: "se",
-    35: "br",
-    36: "kr",
-    37: "rb",
-    38: "sr",
-    39: "y",
-    40: "zr",
-    41: "nb",
-    42: "mo",
-    43: "tc",
-    44: "ru",
-    45: "rh",
-    46: "pd",
-    47: "ag",
-    48: "cd",
-    49: "in",
-    50: "sn",
-    51: "sb",
-    52: "te",
-    53: "i",
-    54: "xe",
-    55: "cs",
-    56: "ba",
-    57: "la",
-    58: "ce",
-    59: "pr",
-    60: "nd",
-    61: "pm",
-    62: "sm",
-    63: "eu",
-    64: "gd",
-    65: "tb",
-    66: "dy",
-    67: "ho",
-    68: "er",
-    69: "tm",
-    70: "yb",
-    71: "lu",
-    72: "hf",
-    73: "ta",
-    74: "w",
-    75: "re",
-    76: "os",
-    77: "ir",
-    78: "pt",
-    79: "au",
-    80: "hg",
-    81: "tl",
-    82: "pb",
-    83: "bi",
-    84: "po",
-    85: "at",
-    86: "rn",
-    87: "fr",
-    88: "ra",
-    89: "ac",
-    90: "th",
-    91: "pa",
-    92: "u",
-    93: "np",
-    94: "pu",
-    95: "am",
-    96: "cm",
-    97: "bk",
-    98: "cf",
-    99: "es",
-    100: "fm",
-    101: "md",
-    102: "no",
-    103: "lr",
-    999: "q",
-    1000: "cav",
-    1000000: "sup",
-    1: "dum",
-}
-
+atomic_Symbol = atomicNumbersReverse
 
 atomicElectronConfigurations = {
     "h": "1s1",
diff --git a/ThermoScreening/thermo/system.py b/ThermoScreening/thermo/system.py
index 9916dda..eed96a8 100644
--- a/ThermoScreening/thermo/system.py
+++ b/ThermoScreening/thermo/system.py
@@ -4,6 +4,8 @@
 
 from beartype.typing import List
 from numpy.exceptions import ComplexWarning
+from PQAnalysis.atomic_system import AtomicSystem as PQAtomicSystem
+from PQAnalysis.core.atom import Atom as PQAtom
 
 from pymatgen.core import Molecule, Structure
 from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
@@ -73,6 +75,13 @@ def _point_group_analyzer(molecule):
         return PointGroupAnalyzer(molecule)
 
 
+def _pq_atomic_system(atoms: List[Atom]) -> PQAtomicSystem:
+    return PQAtomicSystem(
+        atoms=[PQAtom(atom.symbol) for atom in atoms],
+        pos=np.asarray([_real_position(atom.position) for atom in atoms], dtype=float),
+    )
+
+
 def linearity(atoms: List[Atom]) -> bool:
     """
     Checks if the system is linear or non-linear.
@@ -325,7 +334,7 @@ def mass(atoms: List[Atom]) -> float:
     float
         The mass of the system in u.
     """
-    return np.sum([atom.mass for atom in atoms])
+    return _pq_atomic_system(atoms).mass
 
 
 def rotational_group_calc(atoms: List[Atom]) -> str:
@@ -350,7 +359,7 @@ def rotational_group_calc(atoms: List[Atom]) -> str:
     return symb
 
 
-def center_of_mass(atoms: List[Atom], mass: float) -> np.ndarray:
+def center_of_mass(atoms: List[Atom], _mass: float) -> np.ndarray:
     """
     Calculates the center of mass of the system.
 
@@ -363,7 +372,7 @@ def center_of_mass(atoms: List[Atom], mass: float) -> np.ndarray:
     np.ndarray
         The center of mass of the system.
     """
-    return np.sum([atom.mass * atom.position for atom in atoms], axis=0) / mass
+    return _pq_atomic_system(atoms).center_of_mass
 
 
 def imaginary_frequencies(vibrational_frequencies: np.ndarray) -> np.ndarray:
diff --git a/pyproject.toml b/pyproject.toml
index 91e1069..760c7a6 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -56,7 +56,7 @@ thermo = "ThermoScreening.cli.thermo:main"
 [tool.pylint.main]
 fail-under = 7.0
 persistent = false
-py-version = "3.10"
+py-version = "3.12"
 
 [tool.pylint.reports]
 score = true
diff --git a/tests/calculator/test_dftbplus.py b/tests/calculator/test_dftbplus.py
index 69e96e3..ba47eab 100644
--- a/tests/calculator/test_dftbplus.py
+++ b/tests/calculator/test_dftbplus.py
@@ -172,9 +172,10 @@ def fake_read(path, format=None):
     assert calls == [("geo_opt.gen", "gen")]
 
 
-def test_hessian_read_reshapes_matrix(monkeypatch, tmp_path):
+def test_hessian_read_uses_square_matrix(monkeypatch, tmp_path):
     monkeypatch.chdir(tmp_path)
     np.savetxt("hessian.out", np.arange(36, dtype=float).reshape(6, 6))
+    monkeypatch.setattr(dftbplus_module, "_pq_hessian_reader", None)
 
     class FakeAtoms:
         def get_global_number_of_atoms(self):
@@ -190,6 +191,39 @@ def get_global_number_of_atoms(self):
     np.testing.assert_array_equal(result.ravel(), np.arange(36, dtype=float))
 
 
+def test_hessian_reader_uses_pqanalysis_when_available(monkeypatch):
+    calls = []
+
+    def fake_reader(filename):
+        calls.append(filename)
+        return np.eye(3)
+
+    monkeypatch.setattr(dftbplus_module, "_pq_hessian_reader", fake_reader)
+
+    np.testing.assert_array_equal(
+        dftbplus_module._read_hessian_matrix("hessian.out"),
+        np.eye(3),
+    )
+    assert calls == ["hessian.out"]
+
+
+def test_hessian_read_rejects_wrong_matrix_size(monkeypatch):
+    class FakeAtoms:
+        def get_global_number_of_atoms(self):
+            return 2
+
+    hessian = Hessian.__new__(Hessian)
+    hessian.atoms = FakeAtoms()
+    monkeypatch.setattr(
+        dftbplus_module,
+        "_read_hessian_matrix",
+        lambda filename: np.zeros((3, 3)),
+    )
+
+    with pytest.raises(ValueError, match="Hessian matrix size"):
+        hessian.read()
+
+
 def test_modes_initialization_runs_steps(monkeypatch):
     calls = []
 
diff --git a/tests/conftest.py b/tests/conftest.py
index 3af51a0..9a98f37 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -1,39 +1,21 @@
-import pytest
-import os
 import shutil
-# (c) Jakob Gamper
-
-
-@pytest.fixture(scope="function")
-def tmpdir():
-
-    tmpdir = "tmpdir"
-
-    if os.path.exists(tmpdir) and os.path.isdir(tmpdir):
-        shutil.rmtree(tmpdir)
-    os.mkdir(tmpdir)
+from pathlib import Path
 
-    os.chdir(tmpdir)
-
-    yield tmpdir
-
-    os.chdir("..")
-    shutil.rmtree(tmpdir)
+import pytest
 
 
 @pytest.fixture(scope="function")
-def test_with_data_dir(example_dir):
+def tmpdir(tmp_path, monkeypatch):
+    monkeypatch.chdir(tmp_path)
+    return str(tmp_path)
 
-    tmpdir = "tmpdir"
 
-    if os.path.exists(tmpdir) and os.path.isdir(tmpdir):
-        shutil.rmtree(tmpdir)
-
-    shutil.copytree(os.path.join("tests/data", example_dir), tmpdir)
-
-    os.chdir(tmpdir)
+@pytest.fixture(scope="function")
+def test_with_data_dir(example_dir, tmp_path, monkeypatch):
+    source = Path(__file__).parent / "data" / example_dir
+    workdir = tmp_path / example_dir
 
-    yield tmpdir
+    shutil.copytree(source, workdir)
+    monkeypatch.chdir(workdir)
 
-    os.chdir("..")
-    shutil.rmtree(tmpdir)
+    return str(workdir)
diff --git a/tests/thermo/test_api.py b/tests/thermo/test_api.py
index 484e294..ff75bce 100644
--- a/tests/thermo/test_api.py
+++ b/tests/thermo/test_api.py
@@ -73,11 +73,10 @@ def test_public_api_rejects_wrong_argument_type(self):
             
    
     @patch(
-        "builtins.open",
-        new_callable=mock_open,
-        read_data="24\n\nO     0.00000003     -0.00000060     -2.14906255      6.45170211\n    O     -0.00000001      0.00000032     -7.56375986      6.45170117\n    C     -3.70943130     -0.00000041     -5.55541507      4.06573254\n    C     -2.50358224      0.00000055     -6.25344004      4.05635404\n    C     -1.28152402      0.00000064     -5.56329596      4.05319169\n    C     -1.28152226     -0.00000048     -4.14952765      4.05319299\n    C     -2.50358591     -0.00000218     -3.45938372      4.05635329\n    C     -3.70942822     -0.00000196     -4.15740715      4.06573042\n    C      0.00000004      0.00000235     -6.33321909      3.57768382\n    C     -0.00000003      0.00000103     -3.37960436      3.57768340\n    C      1.28152286      0.00000564     -4.14952745      4.05319251\n    C      1.28152346      0.00000677     -5.56329581      4.05319217\n    C      2.50358346      0.00001200     -6.25343979      4.05635377\n    H      2.48494730      0.00001319     -7.34030311      0.89492446\n    C      3.70943026      0.00001569     -5.55541540      4.06573185\n    C      3.70942928      0.00001410     -4.15740743      4.06573112\n    C      2.50358473      0.00000924     -3.45938351      4.05635352\n    H     -4.65151464      0.00000015     -6.09787561      0.91510467\n    H     -2.48494736      0.00000157     -7.34030301      0.89492435\n    H     -2.48494723     -0.00000334     -2.37251966      0.89492456\n    H     -4.65151524     -0.00000279     -3.61494721      0.91510621\n    H      4.65151483      0.00001973     -6.09787570      0.91510518\n    H      4.65151507      0.00001672     -3.61494731      0.91510569\n    H      2.48494732      0.00000821     -2.37251974      0.89492448",
+        "pathlib.Path.read_text",
+        return_value="24\n\nO     0.00000003     -0.00000060     -2.14906255      6.45170211\n    O     -0.00000001      0.00000032     -7.56375986      6.45170117\n    C     -3.70943130     -0.00000041     -5.55541507      4.06573254\n    C     -2.50358224      0.00000055     -6.25344004      4.05635404\n    C     -1.28152402      0.00000064     -5.56329596      4.05319169\n    C     -1.28152226     -0.00000048     -4.14952765      4.05319299\n    C     -2.50358591     -0.00000218     -3.45938372      4.05635329\n    C     -3.70942822     -0.00000196     -4.15740715      4.06573042\n    C      0.00000004      0.00000235     -6.33321909      3.57768382\n    C     -0.00000003      0.00000103     -3.37960436      3.57768340\n    C      1.28152286      0.00000564     -4.14952745      4.05319251\n    C      1.28152346      0.00000677     -5.56329581      4.05319217\n    C      2.50358346      0.00001200     -6.25343979      4.05635377\n    H      2.48494730      0.00001319     -7.34030311      0.89492446\n    C      3.70943026      0.00001569     -5.55541540      4.06573185\n    C      3.70942928      0.00001410     -4.15740743      4.06573112\n    C      2.50358473      0.00000924     -3.45938351      4.05635352\n    H     -4.65151464      0.00000015     -6.09787561      0.91510467\n    H     -2.48494736      0.00000157     -7.34030301      0.89492435\n    H     -2.48494723     -0.00000334     -2.37251966      0.89492456\n    H     -4.65151524     -0.00000279     -3.61494721      0.91510621\n    H      4.65151483      0.00001973     -6.09787570      0.91510518\n    H      4.65151507      0.00001672     -3.61494731      0.91510569\n    H      2.48494732      0.00000821     -2.37251974      0.89492448",
     )
-    def test_read_coord(self,mock_open):
+    def test_read_coord(self, read_text_mock):
         
         coord = np.array([[ 0.00000003  ,   -0.00000060  ,   -2.14906255],
                             [-0.00000001 ,    0.00000032  ,   -7.56375986],
@@ -144,16 +143,15 @@ def test_read_coord(self,mock_open):
         assert pbc == False
 
     @patch(
-        "builtins.open",
-        new_callable=mock_open,
-        read_data=(
+        "pathlib.Path.read_text",
+        return_value=(
             "2 10.0 11.0 12.0 90.0 90.0 90.0\n"
             "\n"
             "Cl 0.0 0.0 0.0\n"
             "Na 1.0 2.0 3.0\n"
         ),
     )
-    def test_read_xyz_keeps_multichar_symbols(self, mock_open):
+    def test_read_xyz_keeps_multichar_symbols(self, read_text_mock):
         data_N, data_atoms, data_xyz, cell, pbc = read_xyz("test.xyz")
 
         assert data_N == 2
@@ -163,22 +161,16 @@ def test_read_xyz_keeps_multichar_symbols(self, mock_open):
         assert pbc is True
 
     @patch(
-        "builtins.open",
-        new_callable=mock_open,
-        read_data=(
+        "pathlib.Path.read_text",
+        return_value=(
             "1 bad 11.0 12.0 90.0 90.0 90.0\n"
             "\n"
             "Cl 0.0 0.0 0.0\n"
         ),
     )
-    def test_read_xyz_ignores_invalid_cell_header(self, mock_open):
-        data_N, data_atoms, data_xyz, cell, pbc = read_xyz("test.xyz")
-
-        assert data_N == 1
-        np.testing.assert_array_equal(data_atoms, np.array(["Cl"], dtype=object))
-        np.testing.assert_allclose(data_xyz, np.array([[0.0, 0.0, 0.0]]))
-        assert cell is None
-        assert pbc is False
+    def test_read_xyz_rejects_invalid_cell_header(self, read_text_mock):
+        with pytest.raises(TSValueError, match="Invalid XYZ coordinate file"):
+            read_xyz("test.xyz")
 
     @patch(
         "builtins.open",
@@ -192,6 +184,16 @@ def test_read_vib_file(self, mock_open):
             data_vib, vibrational_frequencies, decimal=8
         )
 
+    @patch("builtins.open", new_callable=mock_open, read_data="1.0\n")
+    def test_read_vib_file_rejects_missing_frequency_column(self, mock_open):
+        with pytest.raises(TSValueError, match="Invalid vibrational frequency line 1"):
+            read_vibrational("test.vib", "dftb+")
+
+    @patch("builtins.open", new_callable=mock_open, read_data="\n")
+    def test_read_vib_file_rejects_empty_file(self, mock_open):
+        with pytest.raises(TSValueError, match="No vibrational frequencies found"):
+            read_vibrational("test.vib", "dftb+")
+
     def test_read_gen(self):
         gen_file = Path(__file__).resolve().parents[1] / "data/thermo/geo_opt.gen"