From aad288029ac06c1f129c12d7d900dc9a025092ec Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonathan=20Sch=C3=B6ps?= <106986430+jonathan-schoeps@users.noreply.github.com> Date: Mon, 18 Nov 2024 11:35:27 +0100 Subject: [PATCH] Bug fix were uhf was always 0 if a fixed charge were given (#80) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * bug fix were uhf was always 0 if a fixed charge were given Signed-off-by: Jonathan Schöps * updated a test for the right uhf Signed-off-by: Jonathan Schöps * get rid of redundant code within set_random_charge Signed-off-by: Jonathan Schöps * using atlist as a variable in misscellaneous Signed-off-by: Jonathan Schöps --------- Signed-off-by: Jonathan Schöps --- .../molecules/generate_molecule.py | 3 +- src/mindlessgen/molecules/miscellaneous.py | 59 ++++++++----------- src/mindlessgen/molecules/refinement.py | 3 +- test/test_generate/test_generate_molecule.py | 31 ++++++++++ 4 files changed, 61 insertions(+), 35 deletions(-) diff --git a/src/mindlessgen/molecules/generate_molecule.py b/src/mindlessgen/molecules/generate_molecule.py index cc6408f..e003ed0 100644 --- 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 index cdf641e..e63e3aa 100644 --- 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 index 4e884c9..a520ebd 100644 --- 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 index e8f105e..8af3623 100644 --- 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