Added support for Turbomole (grimme-lab#101)

* Interface for the qm program turbomole Signed-off-by: Jonathan Schöps <s6jtscho@uni-bonn.de> * support for turbomole Signed-off-by: Jonathan Schöps <s6jtscho@uni-bonn.de> * corrected some typos Signed-off-by: Jonathan Schöps <s6jtscho@uni-bonn.de> * changed the PARNODES command Signed-off-by: Jonathan Schöps <s6jtscho@uni-bonn.de> * coord handling Signed-off-by: Jonathan Schöps <s6jtscho@uni-bonn.de> * requestet changes are build in Signed-off-by: Jonathan Schöps <s6jtscho@uni-bonn.de> * Implementation of the requested changes and the test of the read_mol_from_coord function Signed-off-by: Jonathan Schöps <s6jtscho@uni-bonn.de> * Implementation of two pathes for turbomole Signed-off-by: Jonathan Schöps <s6jtscho@uni-bonn.de> * Ignore HOMO-LUMO gap not implemented as refine engine for Orca and Trubomole Signed-off-by: Jonathan Schöps <s6jtscho@uni-bonn.de> * improve TURBOMOLE configuration Signed-off-by: Marcel Müller <marcel.mueller@thch.uni-bonn.de> * correct linting errors Signed-off-by: Marcel Mueller <marcel.mueller@thch.uni-bonn.de> * exclude more recent interfaces from coverage check Signed-off-by: Marcel Müller <marcel.mueller@thch.uni-bonn.de> * added a test for the 'write_coord_to_file' function Signed-off-by: Jonathan Schöps <s6jtscho@uni-bonn.de> * sync the lokal branch with the remote branch Signed-off-by: Jonathan Schöps <s6jtscho@uni-bonn.de> --------- Signed-off-by: Jonathan Schöps <s6jtscho@uni-bonn.de> Signed-off-by: Marcel Müller <marcel.mueller@thch.uni-bonn.de> Signed-off-by: Marcel Mueller <marcel.mueller@thch.uni-bonn.de> Co-authored-by: Marcel Müller <marcel.mueller@thch.uni-bonn.de>
jonathan-schoeps · Jan 17, 2025 · 89a3d58 · 89a3d58
1 parent a2684a5
commit 89a3d58
Show file tree

Hide file tree

Showing 11 changed files with 687 additions and 15 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -45,6 +45,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - support for `python-3.13`
 - option to set a fixed molecular charge, while ensuring `uhf = 0`
 - `element_composition` and `forbidden_elements` can now be directly set to a `dict` or `list`, respectively, via API access
+- support for TURBOMOLE as QM engine.
 
 ### Breaking Changes
 - Removal of the `dist_threshold` flag and in the `-toml` file.

diff --git a/mindlessgen.toml b/mindlessgen.toml
@@ -51,7 +51,7 @@ molecular_charge = "none"
 [refine]
 # > Maximum number of fragment optimization cycles. Options: <int>
 max_frag_cycles = 10
-# > Quantum Mechanics (QM) engine to use. Options: 'xtb', 'orca'
+# > Quantum Mechanics (QM) engine to use. Options: 'xtb', 'orca'. 'turbomole'
 engine = "xtb"
 # > HOMO-LUMO gap threshold applied at the end of the refinement step
 hlgap = 0.5
@@ -60,7 +60,7 @@ hlgap = 0.5
 debug = false
 
 [postprocess]
-# > Engine for the post-processing part. Options: 'xtb', 'orca'
+# > Engine for the post-processing part. Options: 'xtb', 'orca', 'turbomole'
 engine = "orca"
 # > Optimize geometry in the post-processing part. If `false`, only a single-point is conducted. Options: <bool>
 optimize = true
@@ -88,3 +88,15 @@ basis = "def2-SVP"
 gridsize = 1
 # > Maximum number of SCF cycles: Options: <int>
 scf_cycles = 100
+
+[turbomole]
+# > Path to the ridft executable. The name `ridft` is automatically searched for. Options: <str | Path>
+ridft_path = "/path/to/ridft"
+# > Path to the jobex executable. The name `jobex` is automatically searched for. Options: <str | Path>
+jobex_path = "/path/to/jobex"
+# > Functional/Method: Options: <str>
+functional = "PBE"
+# > Basis set: Options: <str>
+basis = "def2-SVP"
+# > Maximum number of SCF cycles: Options: <int>
+scf_cycles = 100
diff --git a/pyproject.toml b/pyproject.toml
@@ -59,7 +59,12 @@ markers = "optional: mark test as optional"
 plugins = ["covdefaults"]
 source = ["./src"]
 # Exclude interfaces to external programs from coverage
-omit = ["src/mindlessgen/qm/xtb.py", "src/mindlessgen/qm/orca.py"]
+omit = [
+    "src/mindlessgen/qm/xtb.py",
+    "src/mindlessgen/qm/orca.py",
+    "src/mindlessgen/qm/tm.py",
+    "src/mindlessgen/qm/gxtb.py",
+]
 
 [tool.coverage.report]
 fail_under = 50
diff --git a/src/mindlessgen/generator/main.py b/src/mindlessgen/generator/main.py
@@ -10,7 +10,18 @@
 import warnings
 
 from ..molecules import generate_random_molecule, Molecule
-from ..qm import XTB, get_xtb_path, QMMethod, ORCA, get_orca_path, GXTB, get_gxtb_path
+from ..qm import (
+    XTB,
+    get_xtb_path,
+    QMMethod,
+    ORCA,
+    get_orca_path,
+    Turbomole,
+    get_ridft_path,
+    get_jobex_path,
+    GXTB,
+    get_gxtb_path,
+)
 from ..molecules import iterative_optimization, postprocess_mol
 from ..prog import ConfigManager
 
@@ -46,6 +57,8 @@ def generator(config: ConfigManager) -> tuple[list[Molecule] | None, int]:
         config,
         get_xtb_path,
         get_orca_path,  # g-xTB cannot be used anyway
+        get_ridft_path,
+        get_jobex_path,
     )
 
     if config.general.postprocess:
@@ -54,6 +67,8 @@ def generator(config: ConfigManager) -> tuple[list[Molecule] | None, int]:
             config,
             get_xtb_path,
             get_orca_path,
+            get_ridft_path,
+            get_jobex_path,
             get_gxtb_path,
         )
     else:
@@ -77,12 +92,12 @@ def generator(config: ConfigManager) -> tuple[list[Molecule] | None, int]:
     for molcount in range(config.general.num_molecules):
         # print a decent header for each molecule iteration
         if config.general.verbosity > 0:
-            print(f"\n{'='*80}")
+            print(f"\n{'=' * 80}")
             print(
-                f"{'='*22} Generating molecule {molcount + 1:<4} of "
-                + f"{config.general.num_molecules:<4} {'='*24}"
+                f"{'=' * 22} Generating molecule {molcount + 1:<4} of "
+                + f"{config.general.num_molecules:<4} {'=' * 24}"
             )
-            print(f"{'='*80}")
+            print(f"{'=' * 80}")
         manager = mp.Manager()
         stop_event = manager.Event()
         cycles = range(config.general.max_cycles)
@@ -270,6 +285,8 @@ def setup_engines(
     cfg: ConfigManager,
     xtb_path_func: Callable,
     orca_path_func: Callable,
+    ridft_path_func: Callable,
+    jobex_path_func: Callable,
     gxtb_path_func: Callable | None = None,
 ):
     """
@@ -291,6 +308,20 @@ def setup_engines(
         except ImportError as e:
             raise ImportError("orca not found.") from e
         return ORCA(path, cfg.orca)
+    elif engine_type == "turbomole":
+        try:
+            jobex_path = jobex_path_func(cfg.turbomole.jobex_path)
+            if not jobex_path:
+                raise ImportError("jobex not found.")
+        except ImportError as e:
+            raise ImportError("jobex not found.") from e
+        try:
+            ridft_path = ridft_path_func(cfg.turbomole.ridft_path)
+            if not ridft_path:
+                raise ImportError("ridft not found.")
+        except ImportError as e:
+            raise ImportError("ridft not found.") from e
+        return Turbomole(jobex_path, ridft_path, cfg.turbomole)
     elif engine_type == "gxtb":
         if gxtb_path_func is None:
             raise ImportError("No callable function for determining the g-xTB path.")

diff --git a/src/mindlessgen/molecules/molecule.py b/src/mindlessgen/molecules/molecule.py
@@ -134,6 +134,11 @@
 PSE_NUMBERS: dict[str, int] = {k.lower(): v for v, k in PSE.items()}
 PSE_SYMBOLS: dict[int, str] = {v: k.lower() for v, k in PSE.items()}
 
+BOHR2AA = (
+    0.529177210544  # taken from https://physics.nist.gov/cgi-bin/cuu/Value?bohrrada0
+)
+AA2BOHR = 1 / BOHR2AA
+
 
 class Molecule:
     """
@@ -264,6 +269,42 @@ def read_mol_from_file(file: str | Path) -> Molecule:
         molecule.name = file_path.stem
         return molecule
 
+    @staticmethod
+    def read_mol_from_coord(file: str | Path) -> Molecule:
+        """
+        Read the XYZ coordinates and the charge of the molecule from a 'coord' file.
+        Thereby, generate a completely new molecule object from scratch.
+
+        Can be called like this:
+            from molecule import Molecule
+            # Call the static method using the class name
+            coord_file = "coord"
+            molecule_instance = Molecule.read_mol_from_coord(coord_file)
+            # Now you can use the molecule_instance as needed
+            print(molecule_instance.name)
+
+        :param file: The 'coord' file to read from.
+        :return: A new instance of Molecule with the read data.
+        """
+        molecule = Molecule()
+        if isinstance(file, str):
+            file_path = Path(file).resolve()
+        elif isinstance(file, Path):
+            file_path = file.resolve()
+        else:
+            raise TypeError("String or Path expected.")
+        molecule.read_xyz_from_coord(file_path)
+        if file_path.with_suffix(".CHRG").exists():
+            molecule.read_charge_from_file(file_path.with_suffix(".CHRG"))
+        else:
+            molecule.charge = 0
+        if file_path.with_suffix(".UHF").exists():
+            molecule.read_uhf_from_file(file_path.with_suffix(".UHF"))
+        else:
+            molecule.uhf = 0
+        molecule.name = file_path.stem
+        return molecule
+
     @property
     def name(self) -> str:
         """
@@ -480,7 +521,7 @@ def get_xyz_str(self) -> str:
         xyz_str += commentline
         for i in range(self.num_atoms):
             xyz_str += (
-                f"{PSE[self.ati[i]+1]:<5} "
+                f"{PSE[self.ati[i] + 1]:<5} "
                 + f"{self.xyz[i, 0]:>12.7f} "
                 + f"{self.xyz[i, 1]:>12.7f} "
                 + f"{self.xyz[i, 2]:>12.7f}\n"
@@ -558,6 +599,118 @@ def read_xyz_from_file(self, filename: str | Path) -> None:
                 self.xyz[i, 2] = float(line[3])
                 self.atlist[self.ati[i]] += 1
 
+    def get_coord_str(self) -> str:
+        """
+        Obtain a string with the full 'coord' file information of the molecule.
+        """
+        coord_str = "$coord\n"
+        for i in range(self.num_atoms):
+            coord_str += (
+                f"{self.xyz[i, 0] / BOHR2AA:>12.7f} "
+                + f"{self.xyz[i, 1] / BOHR2AA:>12.7f} "
+                + f"{self.xyz[i, 2] / BOHR2AA:>12.7f} "
+                + f"{PSE[self.ati[i] + 1]:>5}\n"
+            )
+        coord_str += "$end\n"
+        return coord_str
+
+    def write_coord_to_file(self, filename: str | Path | None = None):
+        """
+        Write the 'coord' file of the molecule to a file.
+
+        The layout of the file is as follows:
+        ```
+        $coord
+        <x1> <y1> <z1> <symbol 1>
+        <x2> <y2> <z2> <symbol 2>
+        ...
+        $end
+        ```
+
+        :param filename: The name of the file to write to.
+        """
+        # raise an error if the number of atoms is not set
+        if self._num_atoms is None:
+            raise ValueError("Number of atoms not set.")
+        if not self._ati.size:
+            raise ValueError("Atomic numbers not set.")
+        if not self._xyz.size:
+            raise ValueError("Atomic coordinates not set.")
+
+        if filename:
+            if not isinstance(filename, Path):
+                filename = Path(filename).resolve()
+        else:
+            filename = Path("mlm_" + self.name + ".coord").resolve()
+
+        with open(filename, "w", encoding="utf8") as f:
+            f.write(self.get_coord_str())
+        # if the charge is set, write it to a '.CHRG' file
+        if self._charge is not None and self._charge != 0:
+            with open(filename.with_suffix(".CHRG"), "w", encoding="utf8") as f:
+                f.write(f"{self.charge}\n")
+        # if the UHF is set, write it to a '.UHF' file
+        if self._uhf is not None and self._uhf > 0:
+            with open(filename.with_suffix(".UHF"), "w", encoding="utf8") as f:
+                f.write(f"{self.uhf}\n")
+
+    def read_xyz_from_coord(self, filename: str | Path) -> None:
+        """
+        Read the XYZ coordinates of the molecule from a 'coord' file.
+
+        The layout of the file is as follows:
+        ```
+        $coord
+         0.00000000000000      0.00000000000000      3.60590687077610     u
+         0.00000000000000      0.00000000000000     -3.60590687077610     u
+        $end
+        ... or ...
+        $coord
+            0.00000000000000      0.00000000000000      2.30704866281983  u
+            0.00000000000000      0.00000000000000     -2.30704866281983  u
+        $redundant
+            number_of_atoms             2
+            degrees_of_freedom          1
+            internal_coordinates        1
+            frozen_coordinates          0
+        # definitions of redundant internals
+        1 k  1.0000000000000 stre   1    2           val=   4.61410
+                1 non zero eigenvalues  of BmBt
+                1           2.000000000    1    0
+                1
+        $end
+
+        :param filename: The name of the file to read from.
+        """
+        with open(filename, encoding="utf8") as f:
+            lines = f.readlines()
+            if lines[0] != "$coord\n":
+                raise ValueError("File does not start with '$coord'.")
+            # number of atoms
+            num_atoms = 0
+            for line in lines:
+                if line.startswith("$coord"):
+                    continue
+                if (
+                    line.startswith("$end")
+                    or line.startswith("$redundant")
+                    or line.startswith("$user-defined")
+                ):
+                    break
+                num_atoms += 1
+            self.num_atoms = num_atoms
+            # read the atomic coordinates
+            self.xyz = np.zeros((self.num_atoms, 3))
+            self.ati = np.zeros(self.num_atoms, dtype=int)
+            self.atlist = np.zeros(103, dtype=int)
+            for i in range(self.num_atoms):
+                line_entries = lines[i + 1].split()
+                self.xyz[i, 0] = float(line_entries[0]) * BOHR2AA
+                self.xyz[i, 1] = float(line_entries[1]) * BOHR2AA
+                self.xyz[i, 2] = float(line_entries[2]) * BOHR2AA
+                self.ati[i] = PSE_NUMBERS[line_entries[3].lower()] - 1
+                self.atlist[self.ati[i]] += 1
+
     def read_charge_from_file(self, filename: str | Path):
         """
         Read the charge of the molecule from a file.

diff --git a/src/mindlessgen/molecules/refinement.py b/src/mindlessgen/molecules/refinement.py
@@ -3,6 +3,7 @@
 to obtain finally a valid molecule.
 """
 
+import warnings
 from pathlib import Path
 import networkx as nx  # type: ignore
 import numpy as np
@@ -151,6 +152,8 @@ def iterative_optimization(
         gap_sufficient = engine.check_gap(
             molecule=rev_mol, threshold=config_refine.hlgap, verbosity=verbosity
         )
+    except NotImplementedError:
+        warnings.warn("HOMO-LUMO gap check not implemented with this engine.")
     except RuntimeError as e:
         raise RuntimeError("HOMO-LUMO gap could not be checked.") from e
     if not gap_sufficient:

diff --git a/src/mindlessgen/prog/__init__.py b/src/mindlessgen/prog/__init__.py
@@ -7,6 +7,7 @@
     GeneralConfig,
     XTBConfig,
     ORCAConfig,
+    TURBOMOLEConfig,
     GXTBConfig,
     GenerateConfig,
     RefineConfig,
@@ -18,6 +19,7 @@
     "GeneralConfig",
     "XTBConfig",
     "ORCAConfig",
+    "TURBOMOLEConfig",
     "GXTBConfig",
     "GenerateConfig",
     "RefineConfig",