diff --git a/src/tcutility/job/__init__.py b/src/tcutility/job/__init__.py
index 69d6ae5f..b989974a 100644
--- a/src/tcutility/job/__init__.py
+++ b/src/tcutility/job/__init__.py
@@ -4,5 +4,6 @@
 from .adf import ADFJob, ADFFragmentJob  # noqa: F401
 from .dftb import DFTBJob  # noqa: F401
 from .crest import CRESTJob, QCGJob  # noqa: F401
+from .xtb import XTBJob  # noqa: F401
 from .orca import ORCAJob  # noqa: F401
 from .nmr import NMRJob  # noqa: F401
diff --git a/src/tcutility/job/crest.py b/src/tcutility/job/crest.py
index 38049d53..098a4f30 100644
--- a/src/tcutility/job/crest.py
+++ b/src/tcutility/job/crest.py
@@ -37,7 +37,6 @@ def _setup_job(self):
             'coords.xyz',
             f'-xnam "{self.xtb_path}"',
             '--noreftopo',
-            '-rthr 1'
             f'-c {self._charge}',
             f'-u {self._spinpol}',
             f'-tnmd {self._temp}',
diff --git a/src/tcutility/job/generic.py b/src/tcutility/job/generic.py
index f260fc3c..1399722b 100644
--- a/src/tcutility/job/generic.py
+++ b/src/tcutility/job/generic.py
@@ -10,6 +10,22 @@
 j = os.path.join
 
 
+def _python_path():
+    '''
+    Sometimes it is necessary to have the Python path as some environments don't have its path.
+    This function attempts to find the Python path and returns it.
+    '''
+    python = sp.run('which python', shell=True, capture_output=True).stdout.decode().strip()
+
+    if python == '' or not os.path.exists(python):
+        python = sp.run('which python3', shell=True, capture_output=True).stdout.decode().strip()
+
+    # we default to the python executable 
+    if python == '' or not os.path.exists(python):
+        python = 'python'
+
+    return python
+
 class Job:
     '''This is the base Job class used to build more advanced classes such as :class:`AMSJob <tcutility.job.ams.AMSJob>` and :class:`ORCAJob <tcutility.job.orca.ORCAJob>`.
     The base class contains an empty :class:`Result <tcutility.results.result.Result>` object that holds the settings. 
@@ -130,7 +146,8 @@ def run(self):
             self.add_postamble(f'rm -r {self.workdir}')
 
         for postscript in self._postscripts:
-            self._postambles.append(f'python {postscript[0]} {" ".join(postscript[1])}')
+            self._postambles.append(f'{_python_path()} {postscript[0]} {" ".join(postscript[1])}')
+            
         # setup the job and check if it was successfull
         setup_success = self._setup_job()
 
@@ -199,8 +216,12 @@ def add_postscript(self, script, *args):
 
     def dependency(self, otherjob: 'Job'):
         '''
-        Set a dependency between this job and another job. This means that this job will run after the other job is finished running succesfully.
+        Set a dependency between this job and otherjob. 
+        This means that this job will run after the other job is finished running succesfully.
         '''
+        if otherjob.can_skip:
+            return
+            
         if hasattr(otherjob, 'slurm_job_id'):
             self.sbatch(dependency=f'afterok:{otherjob.slurm_job_id}')
             self.sbatch(kill_on_invalid_dep='Yes')
diff --git a/src/tcutility/job/postscripts/write_converged_geoms.py b/src/tcutility/job/postscripts/write_converged_geoms.py
index c9fcf3ce..f1ce8d1e 100644
--- a/src/tcutility/job/postscripts/write_converged_geoms.py
+++ b/src/tcutility/job/postscripts/write_converged_geoms.py
@@ -24,3 +24,9 @@
 			amv.write('\n')
 
 			molecule.save(mol, f'converged_mols/{i}.xyz', comment=f'E = {energy * constants.HA2KCALMOL:.6f} kcal/mol')
+
+			if energy == max(energies):
+				molecule.save(mol, 'converged_mols/highest_energy.xyz', comment=f'E = {energy * constants.HA2KCALMOL:.6f} kcal/mol')
+
+			if energy == min(energies):
+				molecule.save(mol, 'converged_mols/lowest_energy.xyz', comment=f'E = {energy * constants.HA2KCALMOL:.6f} kcal/mol')
diff --git a/src/tcutility/job/xtb.py b/src/tcutility/job/xtb.py
new file mode 100644
index 00000000..1f7d8d43
--- /dev/null
+++ b/src/tcutility/job/xtb.py
@@ -0,0 +1,158 @@
+from tcutility.job.generic import Job
+from tcutility import spell_check
+import os
+from typing import Union
+
+
+j = os.path.join
+
+
+class XTBJob(Job):
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+        self.xtb_path = 'xtb'
+        self._options = ['coords.xyz']
+        self._scan = {}
+
+    def _setup_job(self):
+        os.makedirs(self.workdir, exist_ok=True)
+
+        self._molecule.write(j(self.workdir, 'coords.xyz'))
+
+        if self._scan:
+            with open(j(self.workdir, 'scan.inp'), 'w+') as scaninp:
+                for key, value in self._scan.items():
+                    scaninp.write(f'${key}\n')
+                    for line in value:
+                        scaninp.write(f'    {line}\n')
+                scaninp.write('$end\n')
+
+        options = ' '.join(self._options)
+        with open(self.runfile_path, 'w+') as runf:
+            runf.write('#!/bin/sh\n\n')  # the shebang is not written by default by ADF
+            runf.write('\n'.join(self._preambles) + '\n\n')
+            runf.write(f'{self.xtb_path} {options}\n')
+            runf.write('\n'.join(self._postambles))
+
+        return True
+
+    def model(self, method: Union[str, int]):
+        '''
+        Set the method used by XTB. This includes GFN0-xTB, GFN1-xTB, GFN2-xTB and GFNFF.
+
+        Args:
+            method: the method to use. Can be specified by its full name, e.g. 'GFN2-xTB', is the same as 'GFN2' or simply 2.
+        '''
+        if isinstance(method, int):
+            if not 0 >= method >= 2:
+                raise ValueError(f'GFN{method}-xTB does not exist. Please choose one of GFN[0, 1, 2]-xTB.')
+            self._options.append(f'--gfn {method}')
+            return
+
+        if method.lower() in ['gfn0', 'gfn0-xtb']:
+            self._options.append('--gfn 0')
+            return
+
+        if method.lower() in ['gfn1', 'gfn1-xtb']:
+            self._options.append('--gfn 1')
+            return
+
+        if method.lower() in ['gfn2', 'gfn2-xtb']:
+            self._options.append('--gfn 2')
+            return
+
+        if method.lower() in ['gfnff', 'ff']:
+            self._options.append('--gfnff')
+            return
+
+        raise ValueError(f'Did not recognize the method {method} for XTB.')
+    
+    def solvent(self, name: str = None, model: str = 'alpb'):
+        '''
+        Model solvation using the ALPB or GBSA model.
+
+        Args:
+            name: the name of the solvent you want to use. Must be ``None``, ``Acetone``, ``Acetonitrile``, ``CHCl3``, ``CS2``, ``DMSO``, ``Ether``, ``H2O``, ``Methanol``, ``THF`` or ``Toluene``.
+            grid_size: the size of the grid used to construct the solvent accessible surface. Must be ``230``, ``974``, ``2030`` or ``5810``.
+        '''
+        spell_check.check(model, ['alpb', 'gbsa'], ignore_case=True)
+        self._options.append(f'--{model} {name}')
+
+    def spin_polarization(self, val: int):
+        '''
+        Set the spin-polarization of the system.
+        '''
+        self._options.append(f'-u {val}')
+
+    def multiplicity(self, val: int):
+        '''
+        Set the multiplicity of the system. If the value is not one the calculation will also be unrestricted.
+        We use the following values:
+        
+        1) singlet
+        2) doublet
+        3) triplet
+        4) ...
+
+        The multiplicity is equal to 2*S+1 for spin-polarization of S.
+        '''
+        self._options.append(f'-u {(val - 1)//2}')
+
+    def charge(self, val: int):
+        '''
+        Set the charge of the system.
+        '''
+        self._options.append(f'-c {val}')
+
+    def vibrations(self, enable: bool = True):
+        self._options.append('--hess')
+
+    def optimization(self, quality: str = 'Normal', calculate_hess: bool = True):
+        '''
+        Do a geometry optimization and calculate the normal modes of the input structure.
+        
+        Args:
+            quality: the convergence criteria of the optimization. 
+                See https://xtb-docs.readthedocs.io/en/latest/optimization.html#optimization-levels.
+            calculate_hess: whether to calculate the Hessian and do a normal mode analysis after optimization.
+        '''
+        spell_check.check(quality, ['crude', 'sloppy', 'loose', 'lax', 'normal', 'tight', 'vtight', 'extreme'], ignore_case=True)
+
+        if calculate_hess:
+            self._options.append(f'--ohess {quality}')
+        else:
+            self._options.append(f'--opt {quality}')
+
+    def PESScan(self, distances: list = [], angles: list = [], dihedrals: list = [], npoints: int = 20, quality: str = 'Normal', mode: str = 'concerted'):
+        '''
+        Set the task of the job to potential energy surface scan (PESScan).
+
+        Args:
+            distances: sequence of tuples or lists containing ``[atom_index1, atom_index2, start, end]``. 
+                Atom indices start at 1. Distances are given in |angstrom|.
+            angles: sequence of tuples or lists containing ``[atom_index1, atom_index2, atom_index3, start, end]``. 
+                Atom indices start at 1. Angles are given in degrees
+            dihedrals: sequence of tuples or lists containing ``[atom_index1, atom_index2, atom_index3, atom_index4, start, end]``. 
+                Atom indices start at 1. Angles are given in degrees
+            npoints: the number of PES points to optimize.
+
+        .. note::
+            Currently we only support generating settings for 1-dimensional PESScans. 
+            We will add support for N-dimensional PESScans later.
+        '''
+        self.optimization(quality=quality, calculate_hess=False)
+        self._options.append('--input scan.inp')
+
+        self._scan.setdefault('scan', [f'mode={mode}'])
+        # self._scan.setdefault('opt', [f'maxcycles=50'])
+        for i, d in enumerate(distances):
+            self._scan['scan'].append(f'distance: {d[0]},{d[1]},{d[2]}; {d[2]}, {d[3]}, {npoints}')
+        for i, a in enumerate(angles, start=i):
+            self._scan['scan'].append(f'angle: {a[0]},{a[1]},{a[2]},{a[3]}; {a[3]}, {a[4]}, {npoints}')
+        for i, a in enumerate(dihedrals, start=i):
+            self._scan['scan'].append(f'dihedral: {a[0]},{a[1]},{a[2]},{a[3]},{a[4]}; {a[4]}, {a[5]}, {npoints}')
+
+    @property
+    def output_mol_path(self):
+        return f'{self.workdir}/xtbopt.xyz'
+
diff --git a/src/tcutility/results/__init__.py b/src/tcutility/results/__init__.py
index 2c8b119e..ea1a6694 100644
--- a/src/tcutility/results/__init__.py
+++ b/src/tcutility/results/__init__.py
@@ -15,7 +15,7 @@
 import pathlib as pl  # noqa: E402
 
 from .. import slurm  # noqa: E402
-from . import adf, ams, cache, dftb, orca  # noqa: E402
+from . import adf, ams, cache, dftb, orca, xtb  # noqa: E402
 
 
 def get_info(calc_dir: str):
@@ -29,6 +29,11 @@ def get_info(calc_dir: str):
     except:  # noqa
         pass
 
+    try:
+        return xtb.get_info(calc_dir)
+    except:  # noqa
+        pass
+
     res = Result()
 
     # if we cannot correctly read the info, we return some generic result object
@@ -93,6 +98,10 @@ def read(calc_dir: Union[str, pl.Path]) -> Result:
     elif ret.engine == "dftb":
         ret.dftb = dftb.get_calc_settings(ret)
         ret.properties = dftb.get_properties(ret)
+    elif ret.engine == "xtb":
+        # ret.xtb = xtb.get_calc_settings(ret)
+        ret.properties = xtb.get_properties(ret)
+        
     elif ret.engine == "orca":
         try:
             ret.orca = orca.get_calc_settings(ret)
@@ -100,7 +109,6 @@ def read(calc_dir: Union[str, pl.Path]) -> Result:
             ret.orca = None
             print('Error reading:', calc_dir)
             raise
-
         try:
             ret.properties = orca.get_properties(ret)
         except:
diff --git a/src/tcutility/results/ams.py b/src/tcutility/results/ams.py
index 242f1789..6b7e7786 100644
--- a/src/tcutility/results/ams.py
+++ b/src/tcutility/results/ams.py
@@ -24,7 +24,7 @@ def get_calc_files(calc_dir: str) -> dict:
     """
     # collect all files in the current directory and subdirectories
     files = []
-    for root, _, files_ in os.walk(calc_dir):
+    for root, _, files_ in os.walk(os.path.abspath(calc_dir)):
         if os.path.split(root)[1].startswith('.'):
             continue
 
@@ -468,6 +468,12 @@ def get_history(calc_dir: str) -> Result:
                     val = reader_ams.read("History", f"{item}({i+1})")
                     ret[item.lower()].append(val)
 
+        if 'converged' not in ret.keys() and ('PESScan', 'HistoryIndices') in reader_ams:
+            ret['converged'] = [False] * ret.number_of_entries
+            for idx in ensure_list(reader_ams.read('PESScan', 'HistoryIndices')):
+                ret['converged'][idx - 1] = True
+
+
     return ret
 
 
diff --git a/src/tcutility/results/dftb.py b/src/tcutility/results/dftb.py
index 413a585e..d79b4e74 100644
--- a/src/tcutility/results/dftb.py
+++ b/src/tcutility/results/dftb.py
@@ -1,5 +1,5 @@
 from tcutility.results import cache, Result
-from tcutility import ensure_list
+from tcutility import ensure_list, constants
 
 
 def get_calc_settings(info: Result) -> Result:
@@ -8,6 +8,10 @@ def get_calc_settings(info: Result) -> Result:
     """
 
     assert info.engine == "dftb", f"This function reads DFTB data, not {info.engine} data"
+
+    if info.input.task == 'PESScan':
+        return Result()
+
     assert "dftb.rkf" in info.files, f'Missing dftb.rkf file, [{", ".join([": ".join(item) for item in info.files.items()])}]'
 
     ret = Result()
@@ -19,18 +23,22 @@ def get_calc_settings(info: Result) -> Result:
 
 
 def get_properties(info: Result) -> Result:
-    """Function to get properties from an ADF calculation.
+    """Function to get properties from an DFTB calculation.
 
     Args:
-        info: Result object containing ADF properties.
+        info: Result object containing DFTB properties.
 
     Returns:
-        :Result object containing properties from the ADF calculation:
+        :Result object containing properties from the DFTB calculation:
 
             - **energy.bond (float)** – bonding energy (|kcal/mol|).
     """
 
     assert info.engine == "dftb", f"This function reads DFTB data, not {info.engine} data"
+
+    if info.input.task == 'PESScan':
+        return Result()
+
     assert "dftb.rkf" in info.files, f'Missing dftb.rkf file, [{", ".join([": ".join(item) for item in info.files.items()])}]'
 
     reader_dftb = cache.get(info.files["dftb.rkf"])
@@ -48,21 +56,33 @@ def get_properties(info: Result) -> Result:
     for i in range(1, number_of_properties + 1):
         subtypes.append(reader_dftb.read("Properties", f"Subtype({i})").strip())
         types.append(reader_dftb.read("Properties", f"Type({i})").strip())
-        values.append(reader_dftb.read("Properties", f"Value({i})"))
+        if types[-1] == 'Energy':
+            values.append(reader_dftb.read("Properties", f"Value({i})") * constants.HA2KCALMOL)
+        else:
+            values.append(reader_dftb.read("Properties", f"Value({i})"))
 
     # then simply add the properties to ret
     for typ, subtyp, value in zip(types, subtypes, values):
         ret[typ.replace(" ", "_")][subtyp] = value
 
+    if ret.energy['DFTB Final Energy']:
+        ret.energy.bond = ret.energy['DFTB Final Energy']
+
     # we also read vibrations
-    ret.vibrations.number_of_modes = reader_dftb.read("Vibrations", "nNormalModes")
-    ret.vibrations.frequencies = ensure_list(reader_dftb.read("Vibrations", "Frequencies[cm-1]"))
-    if ("Vibrations", "Intensities[km/mol]") in reader_dftb:
-        ret.vibrations.intensities = ensure_list(reader_dftb.read("Vibrations", "Intensities[km/mol]"))
-    ret.vibrations.number_of_imag_modes = len([freq for freq in ret.vibrations.frequencies if freq < 0])
-    ret.vibrations.character = "minimum" if ret.vibrations.number_of_imag_modes == 0 else "transitionstate"
-    ret.vibrations.modes = []
-    for i in range(ret.vibrations.number_of_modes):
-        ret.vibrations.modes.append(reader_dftb.read("Vibrations", f"NoWeightNormalMode({i+1})"))
+    if ('Vibrations', 'nNormalModes') in reader_dftb:
+        ret.vibrations.number_of_modes = reader_dftb.read("Vibrations", "nNormalModes")
+        ret.vibrations.frequencies = ensure_list(reader_dftb.read("Vibrations", "Frequencies[cm-1]"))
+        if ("Vibrations", "Intensities[km/mol]") in reader_dftb:
+            ret.vibrations.intensities = ensure_list(reader_dftb.read("Vibrations", "Intensities[km/mol]"))
+        ret.vibrations.number_of_imag_modes = len([freq for freq in ret.vibrations.frequencies if freq < 0])
+        ret.vibrations.character = "minimum" if ret.vibrations.number_of_imag_modes == 0 else "transitionstate"
+        ret.vibrations.modes = []
+        for i in range(ret.vibrations.number_of_modes):
+            ret.vibrations.modes.append(reader_dftb.read("Vibrations", f"NoWeightNormalMode({i+1})"))
+
+    if ("Thermodynamics", "Gibbs free Energy") in reader_dftb:
+        ret.energy.gibbs = reader_dftb.read("Thermodynamics", "Gibbs free Energy") * constants.HA2KCALMOL
+        ret.energy.enthalpy = reader_dftb.read("Thermodynamics", "Enthalpy") * constants.HA2KCALMOL
+        ret.energy.nuclear_internal = reader_dftb.read("Thermodynamics", "Internal Energy total") * constants.HA2KCALMOL
 
     return ret
diff --git a/src/tcutility/results/xtb.py b/src/tcutility/results/xtb.py
new file mode 100644
index 00000000..b882c93d
--- /dev/null
+++ b/src/tcutility/results/xtb.py
@@ -0,0 +1,440 @@
+from tcutility.results import Result
+from tcutility import constants, molecule
+import os
+import numpy as np
+
+
+j = os.path.join
+
+
+def get_calc_files(calc_dir: str) -> Result:
+    """Function that returns files relevant to AMS calculations stored in ``calc_dir``.
+
+    Args:
+        calc_dir: path pointing to the desired calculation
+
+    Returns:
+        Result object containing filenames and paths
+    """
+    # collect all files in the current directory and subdirectories
+    files = []
+    for root, _, files_ in os.walk(calc_dir):
+        files.extend([j(root, file) for file in files_])
+
+    # parse the filenames
+    ret = Result()
+    ret.extra = []
+    ret.root = os.path.abspath(calc_dir)
+    for file in files:
+        if file.endswith('.out') and not file.endswith('g98.out'):
+            ret.out = os.path.abspath(file)
+            continue
+        if file.endswith('xtbopt.xyz'):
+            ret.opt_out = os.path.abspath(file)
+            continue
+        if file.endswith('xtbopt.log'):
+            ret.opt_history = os.path.abspath(file)
+            continue
+        if file.endswith('xtbscan.log'):
+            ret.scan_out = os.path.abspath(file)
+            continue
+        if file.endswith('vibspectrum'):
+            ret.vibspectrum = os.path.abspath(file)
+            continue
+        if file.endswith('hessian'):
+            ret.hessian = os.path.abspath(file)
+            continue
+
+        ret.extra.append(os.path.abspath(file))
+
+    return ret
+
+
+def get_version(info: Result) -> Result:
+    """Function to get the ORCA version used in the calculation.
+
+    Args:
+        info: Result object containing ORCA calculation settings.
+
+    Returns:
+        :Result object containing results about the ORCA version:
+
+            - **full (str)** – the full version string as written by ORCA.
+            - **major (str)** – major ORCA version.
+            - **minor (str)** – minor ORCA version.
+            - **micro (str)** – micro ORCA version.
+    """
+    ret = Result()
+    with open(info.files.out) as out:
+        for line in out.readlines():
+            line = line.strip()
+            if "* xtb version" not in line:
+                continue
+            version = line.split()[3]
+            ret.full = version
+            ret.major = version.split(".")[0]
+            ret.minor = version.split(".")[1]
+            ret.micro = version.split(".")[2]
+            return ret
+
+
+def get_input(info: Result) -> Result:
+    """Function that parses the input file for this ORCA calculation.
+
+    Args:
+        info: Result object containing ORCA calculation settings.
+
+    Returns:
+        :Result object containing information about the calculation input:
+
+            - **main (list[str])** - the main inputs for the calculation. These are the lines that start with a "!".
+            - **sections (Result)** - extra settings added to the calculation. These are the lines that start with a "%" and optionally end with "END" clause.
+            - **system (Result)** - settings related to the molecular system. This includes charge, multiplicity and the coordinates.
+            - **task (str)** - the task that was performed by the calculation, e.g. "SinglePoint", "TransitionStateSearch".
+    """
+    ret = Result()
+    with open(info.files.out) as out:
+        for line in out.readlines():
+            line = line.strip()
+            if line.startswith('program call'):
+                call = line.split(':')[1].strip().split()
+            if line.startswith('coordinate file'):
+                ret.coord_file = line.split(':')[1].strip()
+
+    ret.call = " ".join(call)
+
+    ### TASK
+    ret.task = 'SinglePoint'
+    for option in ['--opt', '--ohess']:
+        # check if we used the option in our call
+        if option not in call:
+            continue
+
+        # if we used it we did a geo-opt
+        ret.task = 'GeometrOptimization'
+        # check if we did a PESScan, it also requires using the --opt option
+        if 'scan_out' in info.files:
+            ret.task = 'PESScan'
+
+        # check if we gave convergence criteria
+        ret.geometry_convergence = 'Normal'
+        option_idx = call.index(option)
+        if not call[option_idx + 1].startswith('-'):
+            ret.geometry_convergence = call[option_idx + 1]
+
+    ### CHARGE
+    ret.charge = 0
+    for option in ['--charge', '-c']:
+        # check if we used the option in our call
+        if option not in call:
+            continue
+        # read the next position in the call
+        option_idx = call.index(option)
+        ret.charge = int(call[option_idx + 1])
+
+    ### SPIN-POLARIZATION
+    ret.spin_polarization = 0
+    for option in ['--uhf', '-u']:
+        # check if we used the option in our call
+        if option not in call:
+            continue
+        # read the next position in the call
+        option_idx = call.index(option)
+        ret.spin_polarization = int(call[option_idx + 1])
+
+    ### SOLVATION
+    ret.solvent = None
+    ret.solvation_model = None
+
+    if '--alpb' in call:
+        option_idx = call.index('--alpb')
+        ret.solvent = call[option_idx + 1]
+        ret.solvation_model = 'ALPB'
+
+    for option in ['-g', '--gbsa']:
+        if option not in call:
+            continue
+        option_idx = call.index(option)
+        ret.solvent = call[option_idx + 1]
+        ret.solvation_model = 'GBSA'
+
+    ### DETAILED INPUT
+    ret.detailed = None
+    if '--input' in call:
+        ret.detailed = Result()
+        option_idx = call.index('--input')
+        inp_file = call[option_idx + 1]
+        for file in info.files.extra:
+            if file.endswith(inp_file):
+                ret.detailed.file = file
+                break
+
+        with open(ret.detailed.file) as infile:
+            content = infile.readlines()
+
+        for line in content:
+            line = line.strip()
+            if line.startswith('$'):
+                curr_section = line[1:]
+                continue
+
+            if ':' in line:
+                key, val = line.split(':')
+                ret.detailed[curr_section].setdefault(key, [])
+                ret.detailed[curr_section][key].append(val.strip())
+
+            elif '=' in line:
+                key, val = line.split('=')
+                ret.detailed[curr_section][key] = val.strip()
+
+    ### MODEL HAMILTONIAN
+    ret.model = 'GFN2-xTB'
+    if '--gfn' in call:
+        option_idx = call.index('--gfn')
+        version = call[option_idx + 1]
+        ret.model = f'GFN{version}-xTB'
+
+    if '--gfnff' in call:
+        ret.model = 'GFNFF'
+
+    return ret
+
+
+def get_calculation_status(info: Result) -> Result:
+    """Function that returns the status of the ORCA calculation described in reader. In case of non-succes it will also give possible reasons for the errors/warnings.
+
+    Args:
+        info: Result object containing ORCA calculation information.
+
+    Returns:
+        :Result object containing information about the calculation status:
+
+            - **fatal (bool)** – True if calculation cannot be analysed correctly, False otherwise
+            - **reasons (list[str])** – list of reasons to explain the status, they can be errors, warnings, etc.
+            - **name (str)** – calculation status written as a string, one of ("SUCCESS", "RUNNING", "UNKNOWN", "SUCCESS(W)", "FAILED")
+            - **code (str)** – calculation status written as a single character, one of ("S", "R", "U", "W" "F")
+    """
+    ret = Result()
+    ret.fatal = True
+    ret.name = None
+    ret.code = None
+    ret.reasons = []
+
+    if "out" not in info.files:
+        ret.reasons.append("Calculation status unknown")
+        ret.name = "UNKNOWN"
+        ret.code = "U"
+        return ret
+
+    with open(info.files.out) as out:
+        lines = out.readlines()
+        if any(['[WARNING] Runtime exception occurred' in line for line in lines]):
+            ret.fatal = True
+            ret.name = "FAILED"
+            ret.code = "F"
+
+            line_index = [i for i, line in enumerate(lines) if '[WARNING] Runtime exception occurred' in line][0]
+            for line in lines[line_index + 1:]:
+                if '##########' in line:
+                    break
+                ret.reasons.append(line.strip())
+            return ret
+
+        if any(["* finished run" in line for line in lines]):
+            ret.fatal = False
+            ret.name = "SUCCESS"
+            ret.code = "S"
+            return ret
+
+    ret.name = "FAILED"
+    ret.code = "F"
+    return ret
+
+
+def get_molecules(info: Result) -> Result:
+    """Function that returns information about the molecules for this calculation.
+
+    Args:
+        info: Result object containing ORCA calculation information.
+
+    Returns:
+        :Result object containing properties from the ORCA calculation:
+
+            - **input (plams.Molecule)** - the input molecule for this calculation.
+            - **output (plams.Molecule)** - the output molecule for this calculation, for example the optimized structure after a geometry optimization.
+            - **number_of_atoms (int)** - the number of atoms in the molecular system.
+    """
+    ret = Result()
+    for file in info.files.extra:
+        if file.endswith(info.input.coord_file):
+            coord_file = file
+            break
+
+    ret.input = molecule.load(coord_file)
+    ret.output = ret.input.copy()
+
+    if 'opt_out' in info.files:
+        ret.output = molecule.load(info.files.opt_out)
+
+    return ret
+
+
+def get_info(calc_dir: str) -> Result:
+    """Function to read useful info about the calculation in ``calc_dir``. Returned information will depend on the type of file that is provided.
+
+    Args:
+        calc_dir: path pointing to the desired calculation.
+
+    Returns:
+        :Result object containing results about the calculation and AMS:
+
+            - **version (Result)** – information about the AMS version used, see :func:`get_version`.
+            - **files (Result)** - paths to files that are important for this calculation.
+            - **input (Result)** - information about the input of this calculation, see :func:`get_input`.
+            - **level (Result)** - information about the level of theory used for this calculation, see :func:`get_level_of_theory`.
+            - **engine (str)** – the engine that was used to perform the calculation, for example 'adf', 'dftb', ...
+            - **status (Result)** – information about calculation status, see :func:`get_calculation_status`.
+            - **molecule (Result)** – information about the input and output molecules and the molecular system in general, see :func:`get_molecules`.
+    """
+    ret = Result()
+
+    ret.engine = "xtb"
+    ret.files = get_calc_files(calc_dir)
+
+    # store the input of the calculation
+    ret.input = get_input(ret)
+
+    # store information about the version of AMS
+    ret.version = get_version(ret)
+
+    # # store the calculation status
+    ret.status = get_calculation_status(ret)
+
+    # # read molecules
+    ret.molecule = get_molecules(ret)
+
+    return ret
+
+
+
+def get_properties(info: Result) -> Result:
+    """Function to get properties from an ORCA calculation.
+
+    Args:
+        info: Result object containing ORCA properties.
+
+    Returns:
+        :Result object containing properties from the ORCA calculation:
+
+            - **energy.bond (float)** – total bonding energy (|kcal/mol|).
+            - **energy.enthalpy (float)** – enthalpy (|kcal/mol|). Only obtained if vibrational modes were calculated.
+            - **energy.entropy (float)** – entropy (|kcal/mol|). Only obtained if vibrational modes were calculated.
+            - **energy.gibbs (float)** – Gibb's free energy (|kcal/mol|). Only obtained if vibrational modes were calculated.
+            - **energy.[method] (float)** - total energy (|kcal/mol|) at a certain level (HF, MP2, CCSD, ...). This is the sum of energy.HF and energy.[method]_corr.
+            - **energy.[method]_corr (float)** - electron correlation energy (|kcal/mol|) at a certain level (HF, MP2, CCSD, ...).
+            - **vibrations.number_of_modes (int)** – number of vibrational modes for this molecule, 3N-5 for non-linear molecules and 3N-6 for linear molecules, where N is the number of atoms.
+            - **vibrations.number_of_imaginary_modes (int)** – number of imaginary vibrational modes for this molecule.
+            - **vibrations.frequencies (list[float])** – vibrational frequencies associated with the vibrational modes, sorted from low to high (|cm-1|).
+            - **vibrations.intensities (list[float])** – vibrational intensities associated with the vibrational modes (|km/mol|). In ORCA, intensities of imaginary modes are set to 0.
+            - **vibrations.modes (list[float])** – list of vibrational modes sorted from low frequency to high frequency.
+            - **vibrations.character (str)** – the PES character of the molecular system. Either "minimum", "transitionstate" or "saddlepoint(n_imag)", for 0, 1, n_imag number of imaginary frequencies.
+            - **t1** - T1 diagnostic for the highest level of correlation chosen. Used to check the validity of the reference wavefunction.
+            - **s2** - expectation value of the :math:`S^2` operator.
+            - **s2_expected** - ideal expectation value of the :math:`S^2` operator.
+            - **spin_contamination** - the amount of spin-contamination observed in this calculation. It is equal to (s2 - s2_expected) / (s2_expected). Ideally this value should be below 0.1.
+    """
+    ret = Result()
+
+    with open(info.files.out) as out:
+        lines = [line.strip() for line in out.readlines()]
+
+    # read some important info about the calculation
+    for line in lines:
+        if "TOTAL ENERGY" in line:
+            ret.energy.bond = float(line.split()[3]) * constants.HA2KCALMOL
+            continue
+
+        if "TOTAL FREE ENERGY" in line:
+            ret.energy.gibbs = float(line.split()[4]) * constants.HA2KCALMOL
+            continue
+
+        if "TOTAL ENTHALPY" in line:
+            ret.energy.enthalpy = float(line.split()[3]) * constants.HA2KCALMOL
+            continue
+
+        if "# frequencies" in line:
+            ret.vibrations.number_of_modes = int(line.split()[3])
+            continue
+
+        if "# imaginary freq." in line:
+            ret.vibrations.number_of_imaginary_modes = int(line.split()[4])
+            continue
+
+        if 'vibrational frequencies' in line:
+            ret.vibrations.frequencies = []
+            continue
+
+        if 'eigval :' in line:
+            ret.vibrations.frequencies.extend([float(freq) for freq in line.split()[2:]])
+
+    if ret.vibrations.frequencies:
+        # get the number of the first vibrational mode, so without the translations and rotations
+        first_mode = len(ret.vibrations.frequencies) - ret.vibrations.number_of_modes + 1
+        # we have to remove the first 5 or 6 frequencies, because they are translation+rotation
+        ret.vibrations.frequencies = ret.vibrations.frequencies[first_mode - 1:]
+
+        # read in the vibspectrum file to get the vibrational intensities:
+        ret.vibrations.intensities = []
+        if info.files.vibspectrum:
+            with open(info.files.vibspectrum) as spec:
+                for line in spec.readlines()[3:-1]:
+                    if line.startswith('$'):
+                        continue
+
+                    if int(line.split()[0]) < first_mode:
+                        continue
+
+                    ret.vibrations.intensities.append(float(line.split()[3]))
+
+        # read in the hessian file to get the normal modes:
+        if info.files.hessian:
+            hessian = []
+            with open(info.files.hessian) as hess:
+                for line in hess.readlines()[1:]:
+                    hessian.extend([float(x) for x in line.split()])
+
+            # number of modes:
+            N = int(np.sqrt(len(hessian)))
+            # square hessian:
+            H = np.reshape(hessian, (N, N))
+
+            # get the atomic masses
+            mol = molecule.load(info.files.opt_out)
+            masses = np.atleast_2d([atom.mass for atom in mol.atoms for _ in range(3)])
+
+            # make the reduced mass matrix
+            mu = np.sqrt(masses * masses.T)
+            # and reduced mass Hessian
+            F = H / mu
+
+            # diagonalize it to get frequencies and modes
+            freqs, modes = np.linalg.eigh(F)
+            # print(modes)
+
+            # from yviewer import viewer
+            # mol.guess_bonds()
+            # viewer.show(mol, molinfo=[{'normalmode': modes[0].reshape(-1, 3)}])
+
+    return ret
+
+
+if __name__ == "__main__":
+    from tcutility import log
+
+    ret = get_info("/Users/yumanhordijk/PhD/ganesh_project/calculations/AlCl3_xtb/opt_1")
+    log.log(ret.files)
+    props = get_properties(ret)
+    print(props)
+
+    # ret = get_info("/Users/yumanhordijk/PhD/ganesh_project/calculations/AlCl3_xtb/scan_1")
+    # print(ret)