Bug fix were uhf was always 0 if a fixed charge were given (#80)

* bug fix were uhf was always 0 if a fixed charge were given Signed-off-by: Jonathan Schöps <[email protected]> * updated a test for the right uhf Signed-off-by: Jonathan Schöps <[email protected]> * get rid of redundant code within set_random_charge Signed-off-by: Jonathan Schöps <[email protected]> * using atlist as a variable in misscellaneous Signed-off-by: Jonathan Schöps <[email protected]> --------- Signed-off-by: Jonathan Schöps <[email protected]>
grimme-lab · Nov 18, 2024 · aad2880 · aad2880
1 parent 8c3d823
commit aad2880
Show file tree

Hide file tree

Showing 4 changed files with 61 additions and 35 deletions.
diff --git a/src/mindlessgen/molecules/generate_molecule.py b/src/mindlessgen/molecules/generate_molecule.py
@@ -18,6 +18,7 @@
     get_lanthanides,
     get_actinides,
     calculate_ligand_electrons,
+    calculate_uhf,
 )
 
 
@@ -52,7 +53,7 @@ def generate_random_molecule(
         )
     if config_generate.molecular_charge is not None:
         mol.charge = config_generate.molecular_charge
-        mol.uhf = 0
+        mol.uhf = calculate_uhf(mol.atlist)
     else:
         mol.charge, mol.uhf = set_random_charge(mol.ati, verbosity)
     mol.set_name_from_formula()

diff --git a/src/mindlessgen/molecules/miscellaneous.py b/src/mindlessgen/molecules/miscellaneous.py
@@ -4,6 +4,8 @@
 
 import numpy as np
 
+from mindlessgen.molecules.molecule import ati_to_atlist
+
 
 def set_random_charge(
     ati: np.ndarray,
@@ -13,9 +15,8 @@ def set_random_charge(
     Set the charge of a molecule so that unpaired electrons are avoided.
     """
     # go through all ati and add its own value + 1 to get the number of protons
-    nel = 0
-    for atom in ati:
-        nel += atom + 1
+    atlist = ati_to_atlist(ati)
+    nel = calculate_protons(atlist)
     if verbosity > 1:
         print(f"Number of protons in molecule: {nel}")
 
@@ -26,36 +27,9 @@ def set_random_charge(
         # -> The ligands are the remaining protons are assumed to be low spin
         uhf = 0
         charge = 0
-        ln_protons = 0
-        ac_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
-            elif atom in get_actinides():
-                if atom < 96:
-                    uhf += atom - 88
-                else:
-                    uhf += 102 - atom
-                ac_protons += (
-                    atom - 3 + 1
-                )  # subtract 3 to get the number of protons in the Ln3+ ion
-        ligand_protons = nel - ln_protons - ac_protons
-        if verbosity > 2:
-            if np.any(np.isin(ati, get_lanthanides())):
-                print(f"Number of protons from Ln^3+ ions: {ln_protons}")
-            if np.any(np.isin(ati, get_actinides())):
-                print(f"Number of protons from Ac^3+ ions: {ac_protons}")
-            print(
-                f"Number of protons from ligands (assuming negative charge): {ligand_protons}"
-            )
-
-        if ligand_protons % 2 == 0:
+        ligand_electrons = calculate_ligand_electrons(atlist, nel)
+        uhf = calculate_uhf(atlist)
+        if ligand_electrons % 2 == 0:
             charge = 0
         else:
             charge = 1
@@ -111,6 +85,25 @@ def calculate_ligand_electrons(natoms: np.ndarray, nel: int) -> int:
     return ligand_electrons
 
 
+def calculate_uhf(atlist: np.ndarray) -> int:
+    """
+    Calculate the number of unpaired electrons in a molecule.
+    """
+    uhf = 0
+    for ati, occurrence in enumerate(atlist):
+        if ati in get_lanthanides():
+            if ati < 64:
+                uhf += (ati - 56) * occurrence
+            else:
+                uhf += (70 - ati) * occurrence
+        elif ati in get_actinides():
+            if ati < 96:
+                uhf += (ati - 88) * occurrence
+            else:
+                uhf += (102 - ati) * occurrence
+    return uhf
+
+
 def get_alkali_metals() -> list[int]:
     """
     Get the atomic numbers of alkali metals.

diff --git a/src/mindlessgen/molecules/refinement.py b/src/mindlessgen/molecules/refinement.py
@@ -13,6 +13,7 @@
     set_random_charge,
     calculate_protons,
     calculate_ligand_electrons,
+    calculate_uhf,
     get_lanthanides,
     get_actinides,
 )
@@ -229,7 +230,7 @@ def detect_fragments(
         # Update the charge of the fragment molecule
         if molecular_charge is not None:
             fragment_molecule.charge = molecular_charge
-            fragment_molecule.uhf = 0
+            fragment_molecule.uhf = calculate_uhf(fragment_molecule.atlist)
         else:
             fragment_molecule.charge, fragment_molecule.uhf = set_random_charge(
                 fragment_molecule.ati, verbosity

diff --git a/test/test_generate/test_generate_molecule.py b/test/test_generate/test_generate_molecule.py
@@ -402,6 +402,7 @@ def test_fixed_charge(
     default_generate_config.molecular_charge = "3"
     default_generate_config.min_num_atoms = 5
     default_generate_config.max_num_atoms = 15
+    default_generate_config.forbidden_elements = "57-71, 89-103"
 
     # Ensure the charge is correctly set
     mol = generate_random_molecule(default_generate_config, verbosity=1)
@@ -496,3 +497,33 @@ def test_fixed_charge_with_lanthanides_2(default_generate_config):
     atom_list = generate_atom_list(default_generate_config, verbosity=1)
     assert atom_list[0] == 0
     assert atom_list[10] % 2 == 0
+
+
+def test_fixed_charge_with_actinides(default_generate_config):
+    """Test the hydrogen correction for a fixed charge"""
+    default_generate_config.molecular_charge = "3"
+    default_generate_config.min_num_atoms = 5
+    default_generate_config.max_num_atoms = 7
+    default_generate_config.element_composition = (
+        "H:1-1, B:1-1, Ne:2-2, P:1-1, Cl:1-1, Es:1-1"
+    )
+    default_generate_config.forbidden_elements = "2-*"
+
+    # Ensure the right ligand correction is applied
+    mol = generate_random_molecule(default_generate_config, verbosity=1)
+    assert mol.uhf == 4
+
+
+def test_fixed_charge_with_lanthanides_and_actinides(default_generate_config):
+    """Test the hydrogen correction for a fixed charge"""
+    default_generate_config.molecular_charge = "3"
+    default_generate_config.min_num_atoms = 5
+    default_generate_config.max_num_atoms = 7
+    default_generate_config.element_composition = (
+        "B:1-1, Ne:2-2, P:1-1, Cl:1-1, Es:1-1, Pr:1-1"
+    )
+    default_generate_config.forbidden_elements = "2-*"
+
+    # Ensure the right ligand correction is applied
+    mol = generate_random_molecule(default_generate_config, verbosity=1)
+    assert mol.uhf == 6