Changed the uhf value when xtb calculates with lanthanides (grimme-la…

…b#60) * Changed the uhf assingement when xtb calculates with lanthanides Signed-off-by: Jonathan Schöps <[email protected]> * The f electrons are now calculated and substractet correctly Signed-off-by: Jonathan Schöps <[email protected]> * Update CHANGELOG.md Co-authored-by: Marcel Mueller <[email protected]> * The f electron removal for lanthanides are now for the singlepoint implemented Signed-off-by: Jonathan Schöps <[email protected]> * check_ligand_uhf is now a function Signed-off-by: Jonathan Schöps <[email protected]> * The ValueError is now inside the check_ligand_uhf function Signed-off-by: Jonathan Schöps <[email protected]> * check_ligand_uhf now returns nothing but raises a ValueError Signed-off-by: Jonathan Schöps <[email protected]> --------- Signed-off-by: Jonathan Schöps <[email protected]> Co-authored-by: Marcel Mueller <[email protected]>
jonathan-schoeps · Oct 10, 2024 · 3ac09a5 · 3ac09a5
1 parent 2d01578
commit 3ac09a5
Show file tree

Hide file tree

Showing 2 changed files with 40 additions and 2 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -24,6 +24,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Breaking Changes
 - Removal of the `dist_threshold` flag and in the `-toml` file.
+- The number of unpaired electrons (`Molecule.uhf`) is now set to 0 if `xtb` is used as `QMMethod` and a lanthanide is within the molecule to match the `f-in-core` approximation.
 
 ## [0.4.0] - 2024-09-19
 ### Changed

diff --git a/src/mindlessgen/qm/xtb.py b/src/mindlessgen/qm/xtb.py
@@ -6,9 +6,11 @@
 from pathlib import Path
 import shutil
 from tempfile import TemporaryDirectory
+import numpy as np
 from ..molecules import Molecule
 from .base import QMMethod
 from ..prog import XTBConfig
+from ..molecules.miscellaneous import get_lanthanides
 
 
 class XTB(QMMethod):
@@ -34,7 +36,13 @@ def optimize(
         """
         Optimize a molecule using xtb.
         """
-
+        if np.any(np.isin(molecule.ati, get_lanthanides())):
+            check_ligand_uhf(molecule.ati, molecule.charge)
+            # Store the original UHF value and set uhf to 0
+            # Justification: xTB does not treat f electrons explicitly.
+            # The remaining openshell system has to be removed.
+            uhf_original = molecule.uhf
+            molecule.uhf = 0
         # Create a unique temporary directory using TemporaryDirectory context manager
         with TemporaryDirectory(prefix="xtb_") as temp_dir:
             temp_path = Path(temp_dir).resolve()
@@ -71,13 +79,22 @@ def optimize(
             # read the optimized molecule
             optimized_molecule = molecule.copy()
             optimized_molecule.read_xyz_from_file(temp_path / "xtbopt.xyz")
-
+            if np.any(np.isin(molecule.ati, get_lanthanides())):
+                # Reset the UHF value to the original value before returning the optimized molecule.
+                optimized_molecule.uhf = uhf_original
             return optimized_molecule
 
     def singlepoint(self, molecule: Molecule, verbosity: int = 1) -> str:
         """
         Perform a single-point calculation using xtb.
         """
+        if np.any(np.isin(molecule.ati, get_lanthanides())):
+            check_ligand_uhf(molecule.ati, molecule.charge)
+            # Store the original UHF value and set uhf to 0
+            # Justification: xTB does not treat f electrons explicitly.
+            # The remaining openshell system has to be removed.
+            uhf_original = molecule.uhf
+            molecule.uhf = 0
 
         # Create a unique temporary directory using TemporaryDirectory context manager
         with TemporaryDirectory(prefix="xtb_") as temp_dir:
@@ -109,6 +126,8 @@ def singlepoint(self, molecule: Molecule, verbosity: int = 1) -> str:
                     f"xtb failed with return code {return_code}:\n{xtb_log_err}"
                 )
 
+            if np.any(np.isin(molecule.ati, get_lanthanides())):
+                molecule.uhf = uhf_original
             return xtb_log_out
 
     def check_gap(
@@ -201,3 +220,21 @@ def get_xtb_path(binary_name: str | Path | None = None) -> Path:
             xtb_path = Path(which_xtb).resolve()
             return xtb_path
     raise ImportError("'xtb' binary could not be found.")
+
+
+def check_ligand_uhf(ati: np.ndarray, charge: int) -> None:
+    """
+    Check if the remaning number of electrons after the f electrons are removed is even.
+    """
+    nel = 0
+    f_electrons = 0
+    for atom in ati:
+        nel += atom + 1
+        if atom in get_lanthanides():
+            f_electrons += atom - 56
+    # Check if the number of the remaning electrons is even
+    test_rhf = nel - f_electrons - charge
+    if not test_rhf % 2 == 0:
+        raise ValueError(
+            "The remaining number of electrons after the f electrons are removed is not even."
+        )