add Ln support by adaptive UHF recognition (#16)

* add Ln support by adaptive UHF recognition

Signed-off-by: Marcel Müller <marcel.mueller@thch.uni-bonn.de>

* write also UHF file for the sake of consistency

Signed-off-by: Marcel Müller <marcel.mueller@thch.uni-bonn.de>

---------

Signed-off-by: Marcel Müller <marcel.mueller@thch.uni-bonn.de>

CHANGELOG.md

@@ -10,6 +10,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Comment line of `.xyz` file contains the total charge and number of unpaired electrons
 - Default ORCA calculation changed from r2SCAN-3c to PBE/def2-SVP
 - `verbosity = 3` always prints full QM output
+- Adapted generation of number of unpaired electrons; thereby, support for Ln's
+- Shifted group / element sorting definitions to miscellaneous
 
 ### Added
 - Optimization via DFT in the post-processing step

src/mindlessgen/molecules/__init__.py

@@ -7,17 +7,19 @@
     generate_random_molecule,
     generate_coordinates,
     generate_atom_list,
+    check_distances,
+)
+from .refinement import iterative_optimization, detect_fragments
+from .postprocess import postprocess_mol
+from .miscellaneous import (
+    set_random_charge,
     get_three_d_metals,
     get_four_d_metals,
     get_five_d_metals,
     get_lanthanides,
     get_alkali_metals,
     get_alkaline_earth_metals,
-    check_distances,
 )
-from .refinement import iterative_optimization, detect_fragments
-from .postprocess import postprocess_mol
-from .miscellaneous import set_random_charge
 
 __all__ = [
     "Molecule",

src/mindlessgen/molecules/generate_molecule.py

@@ -6,7 +6,15 @@
 import numpy as np
 from ..prog import GenerateConfig
 from .molecule import Molecule
-from .miscellaneous import set_random_charge
+from .miscellaneous import (
+    set_random_charge,
+    get_alkali_metals,
+    get_alkaline_earth_metals,
+    get_three_d_metals,
+    get_four_d_metals,
+    get_five_d_metals,
+    get_lanthanides,
+)
 
 
 def generate_random_molecule(
@@ -29,8 +37,7 @@ def generate_random_molecule(
         inc_scaling_factor=config_generate.increase_scaling_factor,
         verbosity=verbosity,
     )
-    mol.charge = set_random_charge(mol.ati, verbosity)
-    mol.uhf = 0
+    mol.charge, mol.uhf = set_random_charge(mol.ati, verbosity)
     mol.set_name_from_formula()
 
     if verbosity > 1:
@@ -359,53 +366,3 @@ def check_distances(xyz: np.ndarray, threshold: float) -> bool:
             if r < threshold:
                 return False
     return True
-
-
-def get_alkali_metals() -> list[int]:
-    """
-    Get the atomic numbers of alkali metals.
-    """
-    alkali = [2, 10, 18, 36, 54]
-    return alkali
-
-
-def get_alkaline_earth_metals() -> list[int]:
-    """
-    Get the atomic numbers of alkaline earth metals.
-    """
-    alkaline = [3, 11, 19, 37, 55]
-    return alkaline
-
-
-def get_three_d_metals() -> list[int]:
-    """
-    Get the atomic numbers of three d metals.
-    """
-    threedmetals = list(range(20, 30))
-
-    return threedmetals
-
-
-def get_four_d_metals() -> list[int]:
-    """
-    Get the atomic numbers of four d metals.
-    """
-
-    fourdmetals = list(range(38, 48))
-    return fourdmetals
-
-
-def get_five_d_metals() -> list[int]:
-    """
-    Get the atomic numbers of five d metals.
-    """
-    fivedmetals = list(range(71, 80))
-    return fivedmetals
-
-
-def get_lanthanides() -> list[int]:
-    """
-    Get the atomic numbers of lanthanides.
-    """
-    lanthanides = list(range(56, 71))
-    return lanthanides

src/mindlessgen/molecules/miscellaneous.py

@@ -5,7 +5,7 @@
 import numpy as np
 
 
-def set_random_charge(ati: np.ndarray, verbosity: int = 1) -> int:
+def set_random_charge(ati: np.ndarray, verbosity: int = 1) -> tuple[int, int]:
     """
     Set the charge of a molecule so that unpaired electrons are avoided.
     """
@@ -15,22 +15,112 @@ def set_random_charge(ati: np.ndarray, verbosity: int = 1) -> int:
         nel += atom + 1
     if verbosity > 1:
         print(f"Number of protons in molecule: {nel}")
-    iseven = False
-    if nel % 2 == 0:
-        iseven = True
-    # if the number of electrons is even, the charge is -2, 0, or 2
-    # if the number of electrons is odd, the charge is -1, 1
-    randint = np.random.rand()
-    if iseven:
-        if randint < 1 / 3:
-            charge = -2
-        elif randint < 2 / 3:
+
+    if np.any(np.isin(ati, get_lanthanides())):
+        ### Special mode for lanthanides
+        # -> always high spin
+        # -> Divide the molecule into Ln3+ ions and negative "ligands"
+        # -> The ligands are the remaining protons are assumed to be low spin
+        uhf = 0
+        charge = 0
+        ln_protons = 0
+        for atom in ati:
+            if atom in get_lanthanides():
+                if atom < 64:
+                    uhf += atom - 56
+                else:
+                    uhf += 70 - atom
+                ln_protons += (
+                    atom - 3 + 1
+                )  # subtract 3 to get the number of protons in the Ln3+ ion
+        ligand_protons = nel - ln_protons
+        if verbosity > 2:
+            print(f"Number of protons from Ln^3+ ions: {ln_protons}")
+            print(
+                f"Number of protons from ligands (assuming negative charge): {ligand_protons}"
+            )
+        if ligand_protons % 2 == 0:
             charge = 0
         else:
-            charge = 2
+            charge = 1
+        return charge, uhf
     else:
-        if randint < 0.5:
-            charge = -1
+        ### Default mode
+        iseven = False
+        if nel % 2 == 0:
+            iseven = True
+        # if the number of electrons is even, the charge is -2, 0, or 2
+        # if the number of electrons is odd, the charge is -1, 1
+        randint = np.random.rand()
+        if iseven:
+            if randint < 1 / 3:
+                charge = -2
+            elif randint < 2 / 3:
+                charge = 0
+            else:
+                charge = 2
         else:
-            charge = 1
-    return charge
+            if randint < 0.5:
+                charge = -1
+            else:
+                charge = 1
+        uhf = 0
+        return charge, uhf
+
+
+def get_alkali_metals() -> list[int]:
+    """
+    Get the atomic numbers of alkali metals.
+    """
+    alkali = [2, 10, 18, 36, 54]
+    return alkali
+
+
+def get_alkaline_earth_metals() -> list[int]:
+    """
+    Get the atomic numbers of alkaline earth metals.
+    """
+    alkaline = [3, 11, 19, 37, 55]
+    return alkaline
+
+
+def get_three_d_metals() -> list[int]:
+    """
+    Get the atomic numbers of three d metals.
+    """
+    threedmetals = list(range(20, 30))
+
+    return threedmetals
+
+
+def get_four_d_metals() -> list[int]:
+    """
+    Get the atomic numbers of four d metals.
+    """
+
+    fourdmetals = list(range(38, 48))
+    return fourdmetals
+
+
+def get_five_d_metals() -> list[int]:
+    """
+    Get the atomic numbers of five d metals.
+    """
+    fivedmetals = list(range(71, 80))
+    return fivedmetals
+
+
+def get_lanthanides() -> list[int]:
+    """
+    Get the atomic numbers of lanthanides.
+    """
+    lanthanides = list(range(56, 71))
+    return lanthanides
+
+
+def get_actinides() -> list[int]:
+    """
+    Get the atomic numbers of actinides.
+    """
+    actinides = list(range(88, 103))
+    return actinides

src/mindlessgen/molecules/molecule.py

@@ -340,6 +340,9 @@ def uhf(self, value: int | float):
         except ValueError as e:
             raise TypeError("Integer expected.") from e
 
+        if value < 0:
+            raise ValueError("Number of unpaired electrons cannot be negative.")
+
         self._uhf = value
 
     @property
@@ -483,6 +486,10 @@ def write_xyz_to_file(self, filename: str | Path | None = None):
         if self._charge is not None:
             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:
+            with open(filename.with_suffix(".UHF"), "w", encoding="utf8") as f:
+                f.write(f"{self.uhf}\n")
 
     def read_xyz_from_file(self, filename: str | Path):
         """

src/mindlessgen/molecules/refinement.py

@@ -159,8 +159,9 @@ def detect_fragments(mol: Molecule, verbosity: int = 1) -> list[Molecule]:
         for atom in fragment_molecule.ati:
             fragment_molecule.atlist[atom] += 1
         # Update the charge of the fragment molecule
-        fragment_molecule.charge = set_random_charge(fragment_molecule.ati, verbosity)
-        fragment_molecule.uhf = 0
+        fragment_molecule.charge, fragment_molecule.uhf = set_random_charge(
+            fragment_molecule.ati, verbosity
+        )
         fragment_molecule.set_name_from_formula()
         if verbosity > 1:
             print(f"Fragment molecule: {fragment_molecule}")

src/mindlessgen/qm/xtb.py

@@ -47,11 +47,14 @@ def optimize(
                 "--opt",
                 "--gfn",
                 "2",
-                "--chrg",
-                str(molecule.charge),
             ]
+            if molecule.charge != 0:
+                arguments += ["--chrg", str(molecule.charge)]
+            if molecule.uhf != 0:
+                arguments += ["--uhf", str(molecule.uhf)]
             if max_cycles is not None:
                 arguments += ["--cycles", str(max_cycles)]
+
             if verbosity > 2:
                 print(f"Running command: {' '.join(arguments)}")
 
@@ -87,10 +90,13 @@ def singlepoint(self, molecule: Molecule, verbosity: int = 1) -> str:
                 "molecule.xyz",
                 "--gfn",
                 "2",
-                "--chrg",
-                str(molecule.charge),
             ]
-            if verbosity > 1:
+            if molecule.charge != 0:
+                arguments += ["--chrg", str(molecule.charge)]
+            if molecule.uhf != 0:
+                arguments += ["--uhf", str(molecule.uhf)]
+
+            if verbosity > 2:
                 print(f"Running command: {' '.join(arguments)}")
 
             xtb_log_out, xtb_log_err, return_code = self._run(

test/test_generate/test_generate_molecule.py

@@ -181,6 +181,7 @@ def test_generate_molecule() -> None:
     """
     # create a ConfigManager object with verbosity set to 0
     config = ConfigManager()
+    config.generate.forbidden_elements = "57-71"
     config.general.verbosity = 0
     mol = generate_random_molecule(config.generate, config.general.verbosity)
 

test/test_main/test_main.py

@@ -8,10 +8,11 @@
 @pytest.mark.optional
 def test_generator():
     config = ConfigManager()
-    config.general.engine = "xtb"
+    config.refine.engine = "xtb"
     config.general.max_cycles = 10000
-    config.general.parallel = 8
-    config.general.verbosity = 0
+    config.general.parallel = 4
+    config.general.verbosity = -1
+    config.general.postprocess = False
 
     molecules, exitcode = generator(config)
     assert exitcode == 0