From 743a91bde4589f633a3ea3d91598ead411bdd3ca Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonathan=20Sch=C3=B6ps?= Date: Fri, 25 Oct 2024 17:05:57 +0200 Subject: [PATCH 01/13] A new function to set a fixed charge value MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Jonathan Schöps --- mindlessgen.toml | 2 + .../molecules/generate_molecule.py | 61 +++++++++++++++++++ src/mindlessgen/prog/config.py | 17 ++++++ 3 files changed, 80 insertions(+) diff --git a/mindlessgen.toml b/mindlessgen.toml index 0dc56bd..1d4fd1f 100644 --- a/mindlessgen.toml +++ b/mindlessgen.toml @@ -42,6 +42,8 @@ element_composition = "C:2-3, H:1-2, O:1-2, N:1-*" # > I.e., if an element is specified in 'element_composition' with an occurrence > 0, it will be added to the molecule anyway. # > Example: forbidden_elements = "18,57-*" forbidden_elements = "57-71, 81-*" +# > Fixed Charge for the molecule. Options or None +fixed_charge = 1 [refine] # > Maximum number of fragment optimization cycles. Options: diff --git a/src/mindlessgen/molecules/generate_molecule.py b/src/mindlessgen/molecules/generate_molecule.py index 96d0e20..e3176ed 100644 --- a/src/mindlessgen/molecules/generate_molecule.py +++ b/src/mindlessgen/molecules/generate_molecule.py @@ -304,6 +304,44 @@ def check_composition(): elif max_count is not None and natoms[elem] > max_count: natoms[elem] = max_count + def fixed_charge_hydrogen_correction(): + """ + Correct the number of hydrogen atoms if the number of electrons is odd. + """ + + odd_nel, num_atoms = check_nel(natoms, cfg) + # If a fixed charge is defined, and the number of electrons is odd, add an hydrogen atom + if odd_nel and 0 in valid_elems: + for elem, count_range in cfg.element_composition.items(): + if elem == 0: + min_count, max_count = count_range + else: + min_count, max_count = cfg.min_num_atoms, cfg.max_num_atoms + if natoms[0] < max_count and num_atoms < cfg.max_num_atoms: + natoms[0] = natoms[0] + 1 + if natoms[0] > max_count: + raise SystemExit("NOT TODAY MY FRIEND") + elif natoms[0] > min_count and num_atoms > cfg.min_num_atoms: + natoms[0] = natoms[0] - 1 + if natoms[0] < min_count: + raise SystemExit("NOT TODAY MY FRIEND") + elif cfg.min_num_atoms == cfg.max_num_atoms: + print("Fixed charge element correction failed. Restarting generation.") + generate_atom_list(cfg, verbosity) + return natoms + + # TODO: Make this function ready and then look for the last case how this is implemented. + def fixed_charge_elem_correction(): + odd_nel, num_atoms = check_nel(natoms, cfg) + + if odd_nel: + temp_count = 0 + print("temp_count:", temp_count) + for atom in natoms: + if atom > 0 and atom % 2 == 0: + natoms[atom] += 1 + break + ### ACTUAL WORKFLOW START ### # Add a random number of atoms of random types add_random(low_lim_default_random, max_lim_default_random, 0, 3) @@ -321,6 +359,10 @@ def check_composition(): check_composition() # Check if the number of atoms is within the defined limits check_min_max_atoms() + if cfg.fixed_charge and 0 in valid_elems: + fixed_charge_hydrogen_correction() + if cfg.fixed_charge: + fixed_charge_elem_correction() ### ACTUAL WORKFLOW END ### return natoms @@ -426,3 +468,22 @@ def check_distances( if r < scale_minimal_distance * sum_radii: return False return True + + +def check_nel(natoms: np.ndarray, cfg: GenerateConfig) -> tuple[bool, int]: + """ + Check if the number of electrons is odd. + """ + nel = 0 + temp_count = 0 + num_atoms = 0 + for atom in natoms: + temp_count += 1 + if atom > 0: + nel += (atom) * temp_count + num_atoms += atom + nel = nel - cfg.fixed_charge + if nel % 2 == 0: + return False, num_atoms + else: + return True, num_atoms diff --git a/src/mindlessgen/prog/config.py b/src/mindlessgen/prog/config.py index 2fd3457..3689f61 100644 --- a/src/mindlessgen/prog/config.py +++ b/src/mindlessgen/prog/config.py @@ -183,6 +183,7 @@ def __init__(self: GenerateConfig) -> None: self._scale_fragment_detection: float = 1.25 self._scale_minimal_distance: float = 0.8 self._contract_coords: bool = False + self._fixed_charge: int | None = None def get_identifier(self) -> str: return "generate" @@ -402,6 +403,22 @@ def contract_coords(self, contract_coords: bool): raise TypeError("Contract coords should be a boolean.") self._contract_coords = contract_coords + @property + def fixed_charge(self): + """ + Get the fixed_charge. + """ + return self._fixed_charge + + @fixed_charge.setter + def fixed_charge(self, fixed_charge: int): + """ + Set the fixed_charge. + """ + if not isinstance(fixed_charge, int): + raise TypeError("Fixed charge should be an integer.") + self._fixed_charge = fixed_charge + class RefineConfig(BaseConfig): """ From 156efe4625e53bc2728e11d48e14048e011cde9a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonathan=20Sch=C3=B6ps?= Date: Wed, 30 Oct 2024 10:04:13 +0100 Subject: [PATCH 02/13] latest changes to the fixed charge method MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Jonathan Schöps --- .../molecules/generate_molecule.py | 77 +++++++++++++++---- test/test_generate/test_generate_molecule.py | 23 ++++++ 2 files changed, 87 insertions(+), 13 deletions(-) diff --git a/src/mindlessgen/molecules/generate_molecule.py b/src/mindlessgen/molecules/generate_molecule.py index e3176ed..50afb62 100644 --- a/src/mindlessgen/molecules/generate_molecule.py +++ b/src/mindlessgen/molecules/generate_molecule.py @@ -310,37 +310,87 @@ def fixed_charge_hydrogen_correction(): """ odd_nel, num_atoms = check_nel(natoms, cfg) + print(f"Odd nel: {odd_nel}, num_atoms: {num_atoms}") # If a fixed charge is defined, and the number of electrons is odd, add an hydrogen atom - if odd_nel and 0 in valid_elems: + if odd_nel: for elem, count_range in cfg.element_composition.items(): if elem == 0: min_count, max_count = count_range + break else: min_count, max_count = cfg.min_num_atoms, cfg.max_num_atoms if natoms[0] < max_count and num_atoms < cfg.max_num_atoms: + print("I got here") natoms[0] = natoms[0] + 1 - if natoms[0] > max_count: - raise SystemExit("NOT TODAY MY FRIEND") elif natoms[0] > min_count and num_atoms > cfg.min_num_atoms: natoms[0] = natoms[0] - 1 - if natoms[0] < min_count: - raise SystemExit("NOT TODAY MY FRIEND") elif cfg.min_num_atoms == cfg.max_num_atoms: - print("Fixed charge element correction failed. Restarting generation.") - generate_atom_list(cfg, verbosity) + fixed_charge_elem_correction() + elif min_count == max_count: + fixed_charge_elem_correction() + elif natoms[0] < max_count and num_atoms == cfg.max_num_atoms: + fixed_charge_elem_correction() + elif natoms[0] > min_count and num_atoms == cfg.min_num_atoms: + fixed_charge_elem_correction() return natoms # TODO: Make this function ready and then look for the last case how this is implemented. def fixed_charge_elem_correction(): + """ + Correct the number of atoms if the number of electrons is odd and hydrogen can not be added. + """ odd_nel, num_atoms = check_nel(natoms, cfg) + # if odd_nel: + # temp_count = 0 + # odd_atom_in_list = False + # for atom in natoms: + # if atom > 0 and temp_count > 0 and temp_count % 2 == 0: + # for elem, count_range in cfg.element_composition.items(): + # if elem == temp_count: + # min_count, max_count = count_range + # break + # else: + # min_count, max_count = cfg.min_num_atoms, cfg.max_num_atoms + # + # if natoms[temp_count] < max_count and num_atoms < cfg.max_num_atoms: + # natoms[temp_count] += 1 + # elif ( + # natoms[temp_count] > min_count and num_atoms > cfg.min_num_atoms + # ): + # natoms[temp_count] -= 1 + # odd_atom_in_list = True + # break + # temp_count += 1 if odd_nel: - temp_count = 0 - print("temp_count:", temp_count) - for atom in natoms: - if atom > 0 and atom % 2 == 0: - natoms[atom] += 1 + print("I got until here") + odd_atoms = [elem for elem in valid_elems if elem % 2 == 0] + if 0 in odd_atoms: + odd_atoms.remove(0) + + no_correction = True + while odd_atoms: + print("odd_atoms: ", odd_atoms) + random_elem = np.random.choice(odd_atoms) + print("random_elem: ", random_elem) + for elem, count_range in cfg.element_composition.items(): + if elem == random_elem: + min_count, max_count = count_range + break + else: + min_count, max_count = cfg.min_num_atoms, cfg.max_num_atoms + if natoms[random_elem] < max_count and num_atoms < cfg.max_num_atoms: + natoms[random_elem] += 1 + no_correction = False + break + elif natoms[random_elem] > min_count and num_atoms > cfg.min_num_atoms: + natoms[random_elem] -= 1 + no_correction = False break + else: + odd_atoms.remove(random_elem) + if no_correction: + raise RuntimeError("No valid odd atoms available for correction.") ### ACTUAL WORKFLOW START ### # Add a random number of atoms of random types @@ -358,10 +408,11 @@ def fixed_charge_elem_correction(): # Check if pre-defined atom type counts are within the defined limits check_composition() # Check if the number of atoms is within the defined limits + print(natoms) check_min_max_atoms() if cfg.fixed_charge and 0 in valid_elems: fixed_charge_hydrogen_correction() - if cfg.fixed_charge: + elif cfg.fixed_charge: fixed_charge_elem_correction() ### ACTUAL WORKFLOW END ### diff --git a/test/test_generate/test_generate_molecule.py b/test/test_generate/test_generate_molecule.py index df96b3a..7f27ca5 100644 --- a/test/test_generate/test_generate_molecule.py +++ b/test/test_generate/test_generate_molecule.py @@ -405,3 +405,26 @@ def test_hydrogen_addition( assert atom_list[0] > 0 else: np.testing.assert_equal(atom_list[0], 0) + + +# Check the atom list extention when a fixed charge is given. +@pytest.mark.parametrize( + "charge, expected_charge", + [ + (0, 0), # Neutral charge + (1, 1), # Positive charge + (-1, -1), # Negative charge + ], + ids=["neutral", "positive", "negative"], +) +def test_generate_atom_list_with_charge( + charge, expected_charge, default_generate_config +): + """Test generate_atom_list when a fixed charge is given.""" + default_generate_config.charge = charge + default_generate_config.min_num_atoms = 5 + default_generate_config.max_num_atoms = 15 + atom_list = generate_atom_list(default_generate_config, verbosity=1) + + # Ensure the charge is correctly set + assert np.sum(atom_list) == expected_charge From c73476bf133569382424333d3339e1199ca65211 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonathan=20Sch=C3=B6ps?= Date: Thu, 31 Oct 2024 13:04:54 +0100 Subject: [PATCH 03/13] Fixed charge and element composition works test have to be added MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Jonathan Schöps --- .../molecules/generate_molecule.py | 102 ++++++------------ src/mindlessgen/molecules/miscellaneous.py | 17 ++- test/test_generate/test_generate_molecule.py | 51 ++++++--- 3 files changed, 85 insertions(+), 85 deletions(-) diff --git a/src/mindlessgen/molecules/generate_molecule.py b/src/mindlessgen/molecules/generate_molecule.py index bac816b..1177e9c 100644 --- a/src/mindlessgen/molecules/generate_molecule.py +++ b/src/mindlessgen/molecules/generate_molecule.py @@ -48,6 +48,7 @@ def generate_random_molecule( scale_minimal_distance=config_generate.scale_minimal_distance, ) mol.charge, mol.uhf = set_random_charge(mol.ati, verbosity) + print(mol.charge, mol.uhf) mol.set_name_from_formula() if verbosity > 1: @@ -308,30 +309,26 @@ def fixed_charge_hydrogen_correction(): """ Correct the number of hydrogen atoms if the number of electrons is odd. """ - odd_nel, num_atoms = check_nel(natoms, cfg) - print(f"Odd nel: {odd_nel}, num_atoms: {num_atoms}") - # If a fixed charge is defined, and the number of electrons is odd, add an hydrogen atom - if odd_nel: - for elem, count_range in cfg.element_composition.items(): - if elem == 0: - min_count, max_count = count_range - break - else: - min_count, max_count = cfg.min_num_atoms, cfg.max_num_atoms - if natoms[0] < max_count and num_atoms < cfg.max_num_atoms: - print("I got here") - natoms[0] = natoms[0] + 1 - elif natoms[0] > min_count and num_atoms > cfg.min_num_atoms: - natoms[0] = natoms[0] - 1 - elif cfg.min_num_atoms == cfg.max_num_atoms: - fixed_charge_elem_correction() - elif min_count == max_count: - fixed_charge_elem_correction() - elif natoms[0] < max_count and num_atoms == cfg.max_num_atoms: - fixed_charge_elem_correction() - elif natoms[0] > min_count and num_atoms == cfg.min_num_atoms: - fixed_charge_elem_correction() + if not odd_nel: + return natoms + print("Odd number of electrons. Trying to correct...") + + min_count, max_count = cfg.element_composition.get(0, (0, cfg.max_num_atoms)) + print(min_count, max_count) + print(natoms[0], max_count, num_atoms, cfg.max_num_atoms) + print(natoms[0] < max_count, num_atoms < cfg.max_num_atoms) + print(natoms[0], min_count, num_atoms, cfg.min_num_atoms) + print(natoms[0] > min_count, num_atoms > cfg.min_num_atoms) + if natoms[0] < max_count and num_atoms < cfg.max_num_atoms: + natoms[0] += 1 + print("Adding hydrogen atom...") + elif natoms[0] > min_count and num_atoms > cfg.min_num_atoms: + natoms[0] -= 1 + print("Removing hydrogen atom...") + else: + fixed_charge_elem_correction() + return natoms # TODO: Make this function ready and then look for the last case how this is implemented. @@ -341,56 +338,24 @@ def fixed_charge_elem_correction(): """ odd_nel, num_atoms = check_nel(natoms, cfg) - # if odd_nel: - # temp_count = 0 - # odd_atom_in_list = False - # for atom in natoms: - # if atom > 0 and temp_count > 0 and temp_count % 2 == 0: - # for elem, count_range in cfg.element_composition.items(): - # if elem == temp_count: - # min_count, max_count = count_range - # break - # else: - # min_count, max_count = cfg.min_num_atoms, cfg.max_num_atoms - # - # if natoms[temp_count] < max_count and num_atoms < cfg.max_num_atoms: - # natoms[temp_count] += 1 - # elif ( - # natoms[temp_count] > min_count and num_atoms > cfg.min_num_atoms - # ): - # natoms[temp_count] -= 1 - # odd_atom_in_list = True - # break - # temp_count += 1 if odd_nel: - print("I got until here") - odd_atoms = [elem for elem in valid_elems if elem % 2 == 0] - if 0 in odd_atoms: - odd_atoms.remove(0) - - no_correction = True - while odd_atoms: - print("odd_atoms: ", odd_atoms) - random_elem = np.random.choice(odd_atoms) - print("random_elem: ", random_elem) - for elem, count_range in cfg.element_composition.items(): - if elem == random_elem: - min_count, max_count = count_range - break - else: - min_count, max_count = cfg.min_num_atoms, cfg.max_num_atoms + print("2 time Odd number of electrons. Trying to correct...") + odd_atoms = [elem for elem in valid_elems if elem % 2 == 0 and elem != 0] + + for random_elem in np.random.permutation(odd_atoms): + min_count, max_count = cfg.element_composition.get( + random_elem, (0, cfg.max_num_atoms) + ) if natoms[random_elem] < max_count and num_atoms < cfg.max_num_atoms: natoms[random_elem] += 1 - no_correction = False - break + print(f"Adding atom type {random_elem}...") + return elif natoms[random_elem] > min_count and num_atoms > cfg.min_num_atoms: natoms[random_elem] -= 1 - no_correction = False - break - else: - odd_atoms.remove(random_elem) - if no_correction: - raise RuntimeError("No valid odd atoms available for correction.") + print(f"Removing atom type {random_elem}...") + return + + generate_random_molecule(cfg, verbosity) ### ACTUAL WORKFLOW START ### # Add a random number of atoms of random types @@ -408,7 +373,6 @@ def fixed_charge_elem_correction(): # Check if pre-defined atom type counts are within the defined limits check_composition() # Check if the number of atoms is within the defined limits - print(natoms) check_min_max_atoms() if cfg.fixed_charge and 0 in valid_elems: fixed_charge_hydrogen_correction() diff --git a/src/mindlessgen/molecules/miscellaneous.py b/src/mindlessgen/molecules/miscellaneous.py index 8ed4424..ce03b0e 100644 --- a/src/mindlessgen/molecules/miscellaneous.py +++ b/src/mindlessgen/molecules/miscellaneous.py @@ -4,8 +4,14 @@ import numpy as np +from ..prog import GenerateConfig -def set_random_charge(ati: np.ndarray, verbosity: int = 1) -> tuple[int, int]: +cfg = GenerateConfig() + + +def set_random_charge( + ati: np.ndarray, fixed_charge=cfg.fixed_charge, verbosity: int = 1 +) -> tuple[int, int]: """ Set the charge of a molecule so that unpaired electrons are avoided. """ @@ -51,12 +57,21 @@ def set_random_charge(ati: np.ndarray, verbosity: int = 1) -> tuple[int, int]: print( f"Number of protons from ligands (assuming negative charge): {ligand_protons}" ) + + if fixed_charge: + charge = fixed_charge + return charge, uhf if ligand_protons % 2 == 0: charge = 0 else: charge = 1 return charge, uhf else: + if fixed_charge: + ### Fixed charge mode + charge = fixed_charge + uhf = 0 + return charge, uhf ### Default mode iseven = False if nel % 2 == 0: diff --git a/test/test_generate/test_generate_molecule.py b/test/test_generate/test_generate_molecule.py index 052a8e3..1c8e586 100644 --- a/test/test_generate/test_generate_molecule.py +++ b/test/test_generate/test_generate_molecule.py @@ -408,23 +408,44 @@ def test_hydrogen_addition( # Check the atom list extention when a fixed charge is given. -@pytest.mark.parametrize( - "charge, expected_charge", - [ - (0, 0), # Neutral charge - (1, 1), # Positive charge - (-1, -1), # Negative charge - ], - ids=["neutral", "positive", "negative"], -) -def test_generate_atom_list_with_charge( - charge, expected_charge, default_generate_config +def test_fixed_charge( + default_generate_config, ): - """Test generate_atom_list when a fixed charge is given.""" - default_generate_config.charge = charge + """Test the right assinged charge when a fixed charge is given""" + default_generate_config.fixed_charge = 3 default_generate_config.min_num_atoms = 5 default_generate_config.max_num_atoms = 15 - atom_list = generate_atom_list(default_generate_config, verbosity=1) # Ensure the charge is correctly set - assert np.sum(atom_list) == expected_charge + mol = generate_random_molecule(default_generate_config, verbosity=1) + assert mol.charge == 3 + + +def test_fixed_charge_and_no_possible_correction( + default_generate_config, +): + """Test the hydrogen correction for a fixed charge""" + default_generate_config.fixed_charge = 3 + default_generate_config.min_num_atoms = 5 + default_generate_config.max_num_atoms = 5 + default_generate_config.forbidden_elements = "1-10, 12-*" + + # Ensure the charge is correctly set + mol = generate_random_molecule(default_generate_config, verbosity=1) + assert mol.charge == 3 + assert mol.uhf == 0 + assert mol.atlist[0] == 0 + assert mol.atlist[10] == 5 + + +def test_fixed_charge_hydrogen_correction(default_generate_config): + """Test the hydrogen correction for a fixed charge""" + default_generate_config.fixed_charge = 3 + default_generate_config.min_num_atoms = 5 + default_generate_config.max_num_atoms = 15 + default_generate_config.fixed_elements = "5:1, 10:2, 15:1, 17:1" + default_generate_config.forbidden_elements = "2-10, 11, 12-*" + + # Ensure the charge is correctly set + atom_list = generate_atom_list(default_generate_config, verbosity=1) + assert np.sum(atom_list[0]) % 2 != 0 From 9a21d236b1159f5db0c5aca5bc133306b0df98e0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonathan=20Sch=C3=B6ps?= Date: Wed, 6 Nov 2024 17:04:25 +0100 Subject: [PATCH 04/13] More modular fixed charge routine MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Jonathan Schöps --- src/mindlessgen/generator/main.py | 6 + .../molecules/generate_molecule.py | 213 ++++++++++++------ src/mindlessgen/molecules/miscellaneous.py | 29 +-- src/mindlessgen/molecules/refinement.py | 18 +- src/mindlessgen/qm/xtb.py | 26 ++- test/test_molecules/test_miscellaneous.py | 5 +- test/test_molecules/test_refinement.py | 7 +- 7 files changed, 212 insertions(+), 92 deletions(-) diff --git a/src/mindlessgen/generator/main.py b/src/mindlessgen/generator/main.py index e0c0df8..8a36a54 100644 --- a/src/mindlessgen/generator/main.py +++ b/src/mindlessgen/generator/main.py @@ -173,6 +173,12 @@ def single_molecule_generator( print(e) stop_event.set() return None + except RuntimeError as e: + if config.general.verbosity > 0: + print(f"Generation failed for cycle {cycle + 1}.") + if config.general.verbosity > 1: + print(e) + return None try: # ____ _ _ _ diff --git a/src/mindlessgen/molecules/generate_molecule.py b/src/mindlessgen/molecules/generate_molecule.py index 1177e9c..58ea701 100644 --- a/src/mindlessgen/molecules/generate_molecule.py +++ b/src/mindlessgen/molecules/generate_molecule.py @@ -15,6 +15,7 @@ get_four_d_metals, get_five_d_metals, get_lanthanides, + get_actinides, ) @@ -47,8 +48,11 @@ def generate_random_molecule( ati=mol.ati, scale_minimal_distance=config_generate.scale_minimal_distance, ) - mol.charge, mol.uhf = set_random_charge(mol.ati, verbosity) - print(mol.charge, mol.uhf) + if config_generate.fixed_charge: + mol.charge = config_generate.fixed_charge + mol.uhf = 0 + else: + mol.charge, mol.uhf = set_random_charge(mol.ati, verbosity) mol.set_name_from_formula() if verbosity > 1: @@ -135,6 +139,12 @@ def generate_atom_list(cfg: GenerateConfig, verbosity: int = 1) -> np.ndarray: ) return natoms + # if cfg.fixed_charge and fixed_elem: + # raise ValueError( + # "Both fixed charge and fixed composition are defined. " + # + "Please only define one of them." + # ) + # Reasoning for the parameters in the following sections: # - The number of the atoms added by default (DefaultRandom + AddOrganicAtoms) # should match the minimum number of atoms in the molecule (if defined). @@ -305,57 +315,7 @@ def check_composition(): elif max_count is not None and natoms[elem] > max_count: natoms[elem] = max_count - def fixed_charge_hydrogen_correction(): - """ - Correct the number of hydrogen atoms if the number of electrons is odd. - """ - odd_nel, num_atoms = check_nel(natoms, cfg) - if not odd_nel: - return natoms - print("Odd number of electrons. Trying to correct...") - - min_count, max_count = cfg.element_composition.get(0, (0, cfg.max_num_atoms)) - print(min_count, max_count) - print(natoms[0], max_count, num_atoms, cfg.max_num_atoms) - print(natoms[0] < max_count, num_atoms < cfg.max_num_atoms) - print(natoms[0], min_count, num_atoms, cfg.min_num_atoms) - print(natoms[0] > min_count, num_atoms > cfg.min_num_atoms) - if natoms[0] < max_count and num_atoms < cfg.max_num_atoms: - natoms[0] += 1 - print("Adding hydrogen atom...") - elif natoms[0] > min_count and num_atoms > cfg.min_num_atoms: - natoms[0] -= 1 - print("Removing hydrogen atom...") - else: - fixed_charge_elem_correction() - - return natoms - # TODO: Make this function ready and then look for the last case how this is implemented. - def fixed_charge_elem_correction(): - """ - Correct the number of atoms if the number of electrons is odd and hydrogen can not be added. - """ - odd_nel, num_atoms = check_nel(natoms, cfg) - - if odd_nel: - print("2 time Odd number of electrons. Trying to correct...") - odd_atoms = [elem for elem in valid_elems if elem % 2 == 0 and elem != 0] - - for random_elem in np.random.permutation(odd_atoms): - min_count, max_count = cfg.element_composition.get( - random_elem, (0, cfg.max_num_atoms) - ) - if natoms[random_elem] < max_count and num_atoms < cfg.max_num_atoms: - natoms[random_elem] += 1 - print(f"Adding atom type {random_elem}...") - return - elif natoms[random_elem] > min_count and num_atoms > cfg.min_num_atoms: - natoms[random_elem] -= 1 - print(f"Removing atom type {random_elem}...") - return - - generate_random_molecule(cfg, verbosity) ### ACTUAL WORKFLOW START ### # Add a random number of atoms of random types @@ -374,10 +334,13 @@ def fixed_charge_elem_correction(): check_composition() # Check if the number of atoms is within the defined limits check_min_max_atoms() - if cfg.fixed_charge and 0 in valid_elems: - fixed_charge_hydrogen_correction() + f_block_elements = True + if cfg.fixed_charge and f_block_elements: + natoms = fixed_charge_f_block_correction(cfg, natoms, valid_elems) + elif cfg.fixed_charge and 0 in valid_elems: + fixed_charge_hydrogen_correction(cfg, natoms, valid_elems) elif cfg.fixed_charge: - fixed_charge_elem_correction() + natoms = fixed_charge_elem_correction(cfg, natoms, valid_elems) ### ACTUAL WORKFLOW END ### return natoms @@ -485,20 +448,142 @@ def check_distances( return True -def check_nel(natoms: np.ndarray, cfg: GenerateConfig) -> tuple[bool, int]: +def check_nel(natoms: np.ndarray, cfg: GenerateConfig) -> tuple[bool, int, int, int]: """ Check if the number of electrons is odd. """ - nel = 0 - temp_count = 0 + # nel = 0 + # temp_count = 0 + # num_atoms = 0 + # for atom in natoms: + # temp_count += 1 + # if atom > 0: + # nel += (atom) * temp_count + # num_atoms += atom + # nel = nel - cfg.fixed_charge + temp_count = -1 + nel = -cfg.fixed_charge + f_electrons = 0 num_atoms = 0 for atom in natoms: temp_count += 1 - if atom > 0: - nel += (atom) * temp_count - num_atoms += atom - nel = nel - cfg.fixed_charge + nel += (atom) * (temp_count + 1) + num_atoms += atom + if ( + atom > 0 + and temp_count in get_lanthanides() + or temp_count in get_actinides() + ): + f_electrons = True + if nel % 2 == 0: - return False, num_atoms + return False, num_atoms, f_electrons, nel else: - return True, num_atoms + return True, num_atoms, f_electrons, nel + + +def fixed_charge_f_block_correction( + cfg: GenerateConfig, natoms: np.ndarray, valid_elems: set[int] +) -> np.ndarray: + """ + Correct the number of electrons if f block elements are included. + """ + odd_nel, num_atoms, f_electrons, nel = check_nel(natoms, cfg) + ligand_protons = 0 + ln_protons = 0 + ac_protons = 0 + temp_count = 0 + for atom in natoms: + temp_count += 1 + if (temp_count - 1) in get_lanthanides(): + ln_protons += atom * ( + temp_count - 3 + ) # subtract 3 to get the number of protons in the Ln3+ ion + elif (temp_count - 1) in get_actinides(): + ac_protons += atom * ( + temp_count - 3 + ) # subtract 3 to get the number of protons in the Ln3+ ion + # charge is already within nel + ligand_protons = nel - ln_protons - ac_protons + if ligand_protons % 2 != 0: + odd_atoms = [ + elem + for elem in valid_elems + if elem % 2 == 0 and elem != 0 # TODO: Look for the charge + ] + for random_elem in np.random.permutation(odd_atoms): + min_count, max_count = cfg.element_composition.get( + random_elem, (0, cfg.max_num_atoms) + ) + if natoms[random_elem] < max_count and num_atoms < cfg.max_num_atoms: + natoms[random_elem] += 1 + print(f"Adding atom type {random_elem}...") + return natoms + elif natoms[random_elem] > min_count and num_atoms > cfg.min_num_atoms: + natoms[random_elem] -= 1 + print(f"Removing atom type {random_elem}...") + return natoms + print(odd_atoms, random_elem, min_count, max_count, num_atoms) + raise RuntimeError("Could not correct the odd number of electrons.") + else: + return natoms + + +def fixed_charge_hydrogen_correction( + cfg: GenerateConfig, + natoms: np.ndarray, + valid_elems: set[int], +) -> np.ndarray: + """ + Correct the number of hydrogen atoms if the number of electrons is odd. + """ + odd_nel, num_atoms, f_electrons, nel = check_nel(natoms, cfg) + + if not odd_nel: + return natoms + print("Odd number of electrons. Trying to correct...") + + min_count, max_count = cfg.element_composition.get( + 0, (0, cfg.max_num_atoms) + ) # TODO: Check if max_num_atoms is set when nothing is defined + if natoms[0] < max_count and num_atoms < cfg.max_num_atoms: + natoms[0] += 1 + print("Adding hydrogen atom...") + elif natoms[0] > min_count and num_atoms > cfg.min_num_atoms: + natoms[0] -= 1 + print("Removing hydrogen atom...") + else: + natoms = fixed_charge_elem_correction(cfg, natoms, valid_elems) + + return natoms + + +def fixed_charge_elem_correction( + cfg: GenerateConfig, natoms: np.ndarray, valid_elems: set[int] +) -> np.ndarray: + """ + Correct the number of atoms if the number of electrons is odd and hydrogen can not be added. + """ + odd_nel, num_atoms, f_electrons, nel = check_nel(natoms, cfg) + if odd_nel: + # All other elements + odd_atoms = [ + elem + for elem in valid_elems + if elem % 2 == 0 and elem != 0 # TODO: Look for the charge + ] + for random_elem in np.random.permutation(odd_atoms): + min_count, max_count = cfg.element_composition.get( + random_elem, (0, cfg.max_num_atoms) + ) + if natoms[random_elem] < max_count and num_atoms < cfg.max_num_atoms: + natoms[random_elem] += 1 + print(f"Adding atom type {random_elem}...") + return natoms + elif natoms[random_elem] > min_count and num_atoms > cfg.min_num_atoms: + natoms[random_elem] -= 1 + print(f"Removing atom type {random_elem}...") + return natoms + + raise RuntimeError("Could not correct the odd number of electrons.") + return natoms diff --git a/src/mindlessgen/molecules/miscellaneous.py b/src/mindlessgen/molecules/miscellaneous.py index ce03b0e..68c88f1 100644 --- a/src/mindlessgen/molecules/miscellaneous.py +++ b/src/mindlessgen/molecules/miscellaneous.py @@ -4,13 +4,10 @@ import numpy as np -from ..prog import GenerateConfig - -cfg = GenerateConfig() - def set_random_charge( - ati: np.ndarray, fixed_charge=cfg.fixed_charge, verbosity: int = 1 + ati: np.ndarray, + verbosity: int = 1, ) -> tuple[int, int]: """ Set the charge of a molecule so that unpaired electrons are avoided. @@ -58,20 +55,12 @@ def set_random_charge( f"Number of protons from ligands (assuming negative charge): {ligand_protons}" ) - if fixed_charge: - charge = fixed_charge - return charge, uhf - if ligand_protons % 2 == 0: + elif ligand_protons % 2 == 0: charge = 0 else: charge = 1 return charge, uhf else: - if fixed_charge: - ### Fixed charge mode - charge = fixed_charge - uhf = 0 - return charge, uhf ### Default mode iseven = False if nel % 2 == 0: @@ -95,6 +84,18 @@ def set_random_charge( return charge, uhf +def calculate_protons(natoms: np.ndarray) -> int: + """ + Calculate the number of protons in a molecule. + """ + protons = 0 + temp_count = 0 + for atom in natoms: + temp_count += 1 + protons += atom * temp_count + return protons + + 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 cc8ad4d..cde4564 100644 --- a/src/mindlessgen/molecules/refinement.py +++ b/src/mindlessgen/molecules/refinement.py @@ -50,6 +50,7 @@ def iterative_optimization( # Detect fragments from the optimized molecule fragmols = detect_fragments( mol=rev_mol, + config_generate=config_generate, vdw_scaling=config_generate.scale_fragment_detection, verbosity=verbosity, ) @@ -109,7 +110,7 @@ def iterative_optimization( f"Element {ati} is overrepresented " + f"in the largest fragment in cycle {cycle + 1}." ) - + # TODO: Check the elektrons if uneven throw the fragment away rev_mol = fragmols[ 0 ] # Set current_mol to the first fragment for the next cycle @@ -136,7 +137,10 @@ def iterative_optimization( def detect_fragments( - mol: Molecule, vdw_scaling: float, verbosity: int = 1 + mol: Molecule, + config_generate: GenerateConfig, + vdw_scaling: float, + verbosity: int = 1, ) -> list[Molecule]: """ Detects fragments in a molecular system based on atom positions and covalent radii. @@ -201,9 +205,13 @@ def detect_fragments( for atom in fragment_molecule.ati: fragment_molecule.atlist[atom] += 1 # Update the charge of the fragment molecule - fragment_molecule.charge, fragment_molecule.uhf = set_random_charge( - fragment_molecule.ati, verbosity - ) + if config_generate.fixed_charge is not None: + fragment_molecule.charge = config_generate.fixed_charge + fragment_molecule.uhf = 0 + else: + 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}") diff --git a/src/mindlessgen/qm/xtb.py b/src/mindlessgen/qm/xtb.py index ecd238f..fd4411d 100644 --- a/src/mindlessgen/qm/xtb.py +++ b/src/mindlessgen/qm/xtb.py @@ -10,7 +10,7 @@ from ..molecules import Molecule from .base import QMMethod from ..prog import XTBConfig -from ..molecules.miscellaneous import get_lanthanides +from ..molecules.miscellaneous import get_lanthanides, get_actinides class XTB(QMMethod): @@ -245,13 +245,27 @@ def check_ligand_uhf(ati: np.ndarray, charge: int) -> None: """ nel = 0 f_electrons = 0 + ligand_protons = 0 + ln_protons = 0 for atom in ati: nel += atom + 1 if atom in get_lanthanides(): - f_electrons += atom - 56 + if atom < 64: + f_electrons += atom - 56 + else: + f_electrons += 70 - atom + ln_protons += atom - 3 + 1 + if atom in get_actinides(): + if atom < 96: + f_electrons += atom - 88 + else: + f_electrons += 102 - atom + ln_protons += atom - 3 + 1 + ligand_protons = nel - ln_protons - charge + print("ligand protons are:", ligand_protons) + # 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." + if not ligand_protons % 2 == 0: + raise SystemExit( + "The number of electrons in the ligand is not even. Please check the input." ) diff --git a/test/test_molecules/test_miscellaneous.py b/test/test_molecules/test_miscellaneous.py index 3c6e328..190c51e 100644 --- a/test/test_molecules/test_miscellaneous.py +++ b/test/test_molecules/test_miscellaneous.py @@ -4,8 +4,11 @@ import pytest import numpy as np +from mindlessgen.prog import GenerateConfig # type: ignore from mindlessgen.molecules.miscellaneous import set_random_charge # type: ignore +cfg = GenerateConfig() + # CAUTION: We use 0-based indexing for atoms and molecules! @pytest.mark.parametrize( @@ -69,6 +72,6 @@ def test_set_random_charge(atom_types, expected_charges, expected_uhf): """ Test the set_random_charge function for both standard and lanthanide modes. """ - charge, uhf = set_random_charge(atom_types) + charge, uhf = set_random_charge(atom_types, verbosity=1) assert charge in expected_charges assert uhf == expected_uhf diff --git a/test/test_molecules/test_refinement.py b/test/test_molecules/test_refinement.py index 1471dfd..d075ffd 100644 --- a/test/test_molecules/test_refinement.py +++ b/test/test_molecules/test_refinement.py @@ -7,7 +7,7 @@ from pathlib import Path import numpy as np import pytest -from mindlessgen.prog import ConfigManager # type: ignore +from mindlessgen.prog import ConfigManager, GenerateConfig # type: ignore from mindlessgen.molecules import detect_fragments # type: ignore from mindlessgen.molecules import Molecule # type: ignore from mindlessgen.molecules import iterative_optimization # type: ignore @@ -107,7 +107,10 @@ def test_detect_fragments_H6O2B2Ne2I1Os1Tl1( Test the detection of fragments in the molecule H2O2B2I1Os1. """ fragmols = detect_fragments( - mol=mol_H6O2B2Ne2I1Os1Tl1, vdw_scaling=1.3333, verbosity=0 + mol=mol_H6O2B2Ne2I1Os1Tl1, + config_generate=GenerateConfig(), + vdw_scaling=1.3333, + verbosity=0, ) assert len(fragmols) == 7 From 078109fe0f40e6988c11b957e4288af3f07fd3b4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonathan=20Sch=C3=B6ps?= Date: Fri, 8 Nov 2024 17:03:49 +0100 Subject: [PATCH 05/13] complete fixed charge routine MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Jonathan Schöps --- mindlessgen.toml | 5 +- src/mindlessgen/cli/cli_parser.py | 14 ++ .../molecules/generate_molecule.py | 228 ++++++++---------- src/mindlessgen/molecules/miscellaneous.py | 33 ++- src/mindlessgen/molecules/refinement.py | 23 +- src/mindlessgen/prog/config.py | 35 ++- src/mindlessgen/qm/xtb.py | 3 +- test/test_generate/test_generate_molecule.py | 64 ++++- 8 files changed, 249 insertions(+), 156 deletions(-) diff --git a/mindlessgen.toml b/mindlessgen.toml index 1d4fd1f..5fc1310 100644 --- a/mindlessgen.toml +++ b/mindlessgen.toml @@ -42,8 +42,9 @@ element_composition = "C:2-3, H:1-2, O:1-2, N:1-*" # > I.e., if an element is specified in 'element_composition' with an occurrence > 0, it will be added to the molecule anyway. # > Example: forbidden_elements = "18,57-*" forbidden_elements = "57-71, 81-*" -# > Fixed Charge for the molecule. Options or None -fixed_charge = 1 +# > Set a charge for the molecule generation. Options: , +set_molecular_charge = false +molecular_charge = 0 [refine] # > Maximum number of fragment optimization cycles. Options: diff --git a/src/mindlessgen/cli/cli_parser.py b/src/mindlessgen/cli/cli_parser.py index 47c8484..f2517db 100644 --- a/src/mindlessgen/cli/cli_parser.py +++ b/src/mindlessgen/cli/cli_parser.py @@ -143,6 +143,18 @@ def cli_parser(argv: Sequence[str] | None = None) -> dict: required=False, help="Contract the coordinates of the molecule after the coordinats generation.", ) + parser.add_argument( + "--set-charge", + type=bool, + required=False, + help="Set the charge of the molecule.", + ) + parser.add_argument( + "--molecular-charge", + type=int, + required=False, + help="Define the charge of the molecule.", + ) ### Refinement arguments ### parser.add_argument( @@ -280,6 +292,8 @@ def cli_parser(argv: Sequence[str] | None = None) -> dict: "scale_fragment_detection": args_dict["scale_fragment_detection"], "scale_minimal_distance": args_dict["scale_minimal_distance"], "contract_coords": args_dict["contract_coords"], + "set_molecular_charge": args_dict["set_charge"], + "molecular_charge": args_dict["molecular_charge"], } # XTB specific arguments rev_args_dict["xtb"] = {"xtb_path": args_dict["xtb_path"]} diff --git a/src/mindlessgen/molecules/generate_molecule.py b/src/mindlessgen/molecules/generate_molecule.py index 58ea701..c9fa746 100644 --- a/src/mindlessgen/molecules/generate_molecule.py +++ b/src/mindlessgen/molecules/generate_molecule.py @@ -9,13 +9,14 @@ from .refinement import get_cov_radii, COV_RADII from .miscellaneous import ( set_random_charge, + calculate_protons, get_alkali_metals, get_alkaline_earth_metals, get_three_d_metals, get_four_d_metals, get_five_d_metals, get_lanthanides, - get_actinides, + calculate_ligand_protons, ) @@ -48,8 +49,8 @@ def generate_random_molecule( ati=mol.ati, scale_minimal_distance=config_generate.scale_minimal_distance, ) - if config_generate.fixed_charge: - mol.charge = config_generate.fixed_charge + if config_generate.set_molecular_charge: + mol.charge = config_generate.molecular_charge mol.uhf = 0 else: mol.charge, mol.uhf = set_random_charge(mol.ati, verbosity) @@ -137,14 +138,25 @@ def generate_atom_list(cfg: GenerateConfig, verbosity: int = 1) -> np.ndarray: + "Setting the number of atoms for the fixed composition. " + f"Returning: \n{natoms}\n" ) + # If the molecular charge is defined, and a fixed element composition is defined, check if the electrons are even. If not raise an error. + if cfg.set_molecular_charge: + protons, num_atoms = calculate_protons(natoms) + nel = protons - cfg.molecular_charge + f_elem, ligand_protons = calculate_ligand_protons(natoms, nel) + if f_elem and ligand_protons % 2 != 0 and fixed_elem: + raise SystemExit( + "Both fixed charge and fixed composition are defined. " + + "Please only define one of them." + + "Or ensure that the fixed composition is closed shell." + ) + elif nel % 2 != 0 and fixed_elem: + raise SystemExit( + "Both fixed charge and fixed composition are defined. " + + "Please only define one of them." + + "Or ensure that the fixed composition is closed shell." + ) return natoms - # if cfg.fixed_charge and fixed_elem: - # raise ValueError( - # "Both fixed charge and fixed composition are defined. " - # + "Please only define one of them." - # ) - # Reasoning for the parameters in the following sections: # - The number of the atoms added by default (DefaultRandom + AddOrganicAtoms) # should match the minimum number of atoms in the molecule (if defined). @@ -315,8 +327,6 @@ def check_composition(): elif max_count is not None and natoms[elem] > max_count: natoms[elem] = max_count - # TODO: Make this function ready and then look for the last case how this is implemented. - ### ACTUAL WORKFLOW START ### # Add a random number of atoms of random types add_random(low_lim_default_random, max_lim_default_random, 0, 3) @@ -334,13 +344,13 @@ def check_composition(): check_composition() # Check if the number of atoms is within the defined limits check_min_max_atoms() - f_block_elements = True - if cfg.fixed_charge and f_block_elements: - natoms = fixed_charge_f_block_correction(cfg, natoms, valid_elems) - elif cfg.fixed_charge and 0 in valid_elems: - fixed_charge_hydrogen_correction(cfg, natoms, valid_elems) - elif cfg.fixed_charge: - natoms = fixed_charge_elem_correction(cfg, natoms, valid_elems) + # If f block elements are included, correct the number of electrons + if cfg.set_molecular_charge: + protons, num_atoms = calculate_protons(natoms) + nel = protons - cfg.molecular_charge + natoms = fixed_charge_f_block_correction( + cfg, natoms, num_atoms, nel, valid_elems, verbosity + ) ### ACTUAL WORKFLOW END ### return natoms @@ -448,142 +458,104 @@ def check_distances( return True -def check_nel(natoms: np.ndarray, cfg: GenerateConfig) -> tuple[bool, int, int, int]: - """ - Check if the number of electrons is odd. - """ - # nel = 0 - # temp_count = 0 - # num_atoms = 0 - # for atom in natoms: - # temp_count += 1 - # if atom > 0: - # nel += (atom) * temp_count - # num_atoms += atom - # nel = nel - cfg.fixed_charge - temp_count = -1 - nel = -cfg.fixed_charge - f_electrons = 0 - num_atoms = 0 - for atom in natoms: - temp_count += 1 - nel += (atom) * (temp_count + 1) - num_atoms += atom - if ( - atom > 0 - and temp_count in get_lanthanides() - or temp_count in get_actinides() - ): - f_electrons = True - - if nel % 2 == 0: - return False, num_atoms, f_electrons, nel - else: - return True, num_atoms, f_electrons, nel - - def fixed_charge_f_block_correction( - cfg: GenerateConfig, natoms: np.ndarray, valid_elems: set[int] + cfg: GenerateConfig, + natoms: np.ndarray, + num_atoms: int, + nel: int, + valid_elems: set[int], + verbosity: int, ) -> np.ndarray: """ Correct the number of electrons if f block elements are included. """ - odd_nel, num_atoms, f_electrons, nel = check_nel(natoms, cfg) - ligand_protons = 0 - ln_protons = 0 - ac_protons = 0 - temp_count = 0 - for atom in natoms: - temp_count += 1 - if (temp_count - 1) in get_lanthanides(): - ln_protons += atom * ( - temp_count - 3 - ) # subtract 3 to get the number of protons in the Ln3+ ion - elif (temp_count - 1) in get_actinides(): - ac_protons += atom * ( - temp_count - 3 - ) # subtract 3 to get the number of protons in the Ln3+ ion - # charge is already within nel - ligand_protons = nel - ln_protons - ac_protons - if ligand_protons % 2 != 0: - odd_atoms = [ - elem - for elem in valid_elems - if elem % 2 == 0 and elem != 0 # TODO: Look for the charge - ] - for random_elem in np.random.permutation(odd_atoms): - min_count, max_count = cfg.element_composition.get( - random_elem, (0, cfg.max_num_atoms) + f_elem, ligand_protons = calculate_ligand_protons(natoms, nel) + # If f block elements are included, cerrect only if the remaning ligand protons are uneven + if f_elem and ligand_protons % 2 != 0: + if 0 in valid_elems: + natoms = fixed_charge_hydrogen_correction( + cfg, natoms, num_atoms, valid_elems, verbosity ) - if natoms[random_elem] < max_count and num_atoms < cfg.max_num_atoms: - natoms[random_elem] += 1 - print(f"Adding atom type {random_elem}...") - return natoms - elif natoms[random_elem] > min_count and num_atoms > cfg.min_num_atoms: - natoms[random_elem] -= 1 - print(f"Removing atom type {random_elem}...") - return natoms - print(odd_atoms, random_elem, min_count, max_count, num_atoms) - raise RuntimeError("Could not correct the odd number of electrons.") - else: - return natoms + return natoms + else: + natoms = fixed_charge_elem_correction( + cfg, natoms, num_atoms, valid_elems, verbosity + ) + return natoms + # If f block elements are not included, correct if the number of electrons is uneven + elif nel % 2 != 0 and f_elem is False: + if 0 in valid_elems: + natoms = fixed_charge_hydrogen_correction( + cfg, natoms, num_atoms, valid_elems, verbosity + ) + return natoms + else: + natoms = fixed_charge_elem_correction( + cfg, natoms, num_atoms, valid_elems, verbosity + ) + return natoms + return natoms def fixed_charge_hydrogen_correction( cfg: GenerateConfig, natoms: np.ndarray, + num_atoms: int, valid_elems: set[int], + verbosity: int, ) -> np.ndarray: """ Correct the number of hydrogen atoms if the number of electrons is odd. """ - odd_nel, num_atoms, f_electrons, nel = check_nel(natoms, cfg) - - if not odd_nel: - return natoms - print("Odd number of electrons. Trying to correct...") - - min_count, max_count = cfg.element_composition.get( - 0, (0, cfg.max_num_atoms) - ) # TODO: Check if max_num_atoms is set when nothing is defined + min_count, max_count = cfg.element_composition.get(0, (0, cfg.max_num_atoms)) + if max_count is None: + max_count = cfg.max_num_atoms + # Check if adding or removing hydrogen atoms is possible if natoms[0] < max_count and num_atoms < cfg.max_num_atoms: natoms[0] += 1 - print("Adding hydrogen atom...") + if verbosity > 1: + print("Adding hydrogen atom...") + return natoms elif natoms[0] > min_count and num_atoms > cfg.min_num_atoms: natoms[0] -= 1 - print("Removing hydrogen atom...") + if verbosity > 1: + print("Removing hydrogen atom...") + return natoms else: - natoms = fixed_charge_elem_correction(cfg, natoms, valid_elems) - - return natoms + natoms = fixed_charge_elem_correction( + cfg, natoms, num_atoms, valid_elems, verbosity + ) + return natoms def fixed_charge_elem_correction( - cfg: GenerateConfig, natoms: np.ndarray, valid_elems: set[int] + cfg: GenerateConfig, + natoms: np.ndarray, + num_atoms: int, + valid_elems: set[int], + verbosity: int, ) -> np.ndarray: """ Correct the number of atoms if the number of electrons is odd and hydrogen can not be added. """ - odd_nel, num_atoms, f_electrons, nel = check_nel(natoms, cfg) - if odd_nel: - # All other elements - odd_atoms = [ - elem - for elem in valid_elems - if elem % 2 == 0 and elem != 0 # TODO: Look for the charge - ] - for random_elem in np.random.permutation(odd_atoms): - min_count, max_count = cfg.element_composition.get( - random_elem, (0, cfg.max_num_atoms) - ) - if natoms[random_elem] < max_count and num_atoms < cfg.max_num_atoms: - natoms[random_elem] += 1 - print(f"Adding atom type {random_elem}...") - return natoms - elif natoms[random_elem] > min_count and num_atoms > cfg.min_num_atoms: - natoms[random_elem] -= 1 - print(f"Removing atom type {random_elem}...") - return natoms - - raise RuntimeError("Could not correct the odd number of electrons.") + # All other elements + odd_atoms = [elem for elem in valid_elems if elem % 2 == 0 and elem != 0] + for random_elem in np.random.permutation(odd_atoms): + min_count, max_count = cfg.element_composition.get( + random_elem, (0, cfg.max_num_atoms) + ) + if max_count is None: + max_count = cfg.max_num_atoms + # Check if adding or removing the random element is possible + if natoms[random_elem] < max_count and num_atoms < cfg.max_num_atoms: + natoms[random_elem] += 1 + if verbosity > 1: + print(f"Adding atom type {random_elem} for charge...") + return natoms + elif natoms[random_elem] > min_count and num_atoms > cfg.min_num_atoms: + natoms[random_elem] -= 1 + if verbosity > 1: + print(f"Removing atom type {random_elem} for charge...") + return natoms + raise RuntimeError("Could not correct the odd number of electrons.") return natoms diff --git a/src/mindlessgen/molecules/miscellaneous.py b/src/mindlessgen/molecules/miscellaneous.py index 68c88f1..8d3747e 100644 --- a/src/mindlessgen/molecules/miscellaneous.py +++ b/src/mindlessgen/molecules/miscellaneous.py @@ -84,16 +84,43 @@ def set_random_charge( return charge, uhf -def calculate_protons(natoms: np.ndarray) -> int: +def calculate_protons(natoms: np.ndarray) -> tuple[int, int]: """ - Calculate the number of protons in a molecule. + Calculate the number of protons in a molecule from the atom list. """ protons = 0 + num_atoms = 0 temp_count = 0 for atom in natoms: temp_count += 1 protons += atom * temp_count - return protons + if atom > 0: + num_atoms += atom + return protons, num_atoms + + +def calculate_ligand_protons(natoms: np.ndarray, nel) -> tuple[bool, int]: + """ + Calculate the number of ligand protons in a molecule if lanthanides or actinides are within the molecule. + """ + ligand_protons = 0 + f_protons = 0 + f_elem = False + temp_count = -1 + for atom in natoms: + temp_count += 1 + if ( + atom > 0 + and temp_count in get_lanthanides() + or atom > 0 + and temp_count in get_actinides() + ): + f_protons += atom * ( + temp_count - 3 + 1 + ) # subtract 3 to get the number of protons in the Ln3+ ion + f_elem = True + ligand_protons = nel - f_protons + return f_elem, ligand_protons def get_alkali_metals() -> list[int]: diff --git a/src/mindlessgen/molecules/refinement.py b/src/mindlessgen/molecules/refinement.py index cde4564..fc829c3 100644 --- a/src/mindlessgen/molecules/refinement.py +++ b/src/mindlessgen/molecules/refinement.py @@ -9,7 +9,11 @@ from ..qm.base import QMMethod from ..prog import GenerateConfig, RefineConfig from .molecule import Molecule -from .miscellaneous import set_random_charge +from .miscellaneous import ( + set_random_charge, + calculate_protons, + calculate_ligand_protons, +) COV_RADII = "pyykko" BOHR2AA = ( @@ -110,7 +114,18 @@ def iterative_optimization( f"Element {ati} is overrepresented " + f"in the largest fragment in cycle {cycle + 1}." ) - # TODO: Check the elektrons if uneven throw the fragment away + if config_generate.set_molecular_charge: + protons, num_atoms = calculate_protons(fragmols[0].atlist) + nel = protons - config_generate.molecular_charge + f_elem, ligand_protons = calculate_ligand_protons(fragmols[0].atlist, nel) + if ligand_protons % 2 != 0: + raise RuntimeError( + f"Number of electrons in the largest fragment in cycle {cycle + 1} is odd." + ) + elif nel % 2 != 0: + raise RuntimeError( + f"Number of electrons in the largest fragment in cycle {cycle + 1} is odd." + ) rev_mol = fragmols[ 0 ] # Set current_mol to the first fragment for the next cycle @@ -205,8 +220,8 @@ def detect_fragments( for atom in fragment_molecule.ati: fragment_molecule.atlist[atom] += 1 # Update the charge of the fragment molecule - if config_generate.fixed_charge is not None: - fragment_molecule.charge = config_generate.fixed_charge + if config_generate.set_molecular_charge: + fragment_molecule.charge = config_generate.molecular_charge fragment_molecule.uhf = 0 else: fragment_molecule.charge, fragment_molecule.uhf = set_random_charge( diff --git a/src/mindlessgen/prog/config.py b/src/mindlessgen/prog/config.py index 904a776..931b5ba 100644 --- a/src/mindlessgen/prog/config.py +++ b/src/mindlessgen/prog/config.py @@ -183,7 +183,8 @@ def __init__(self: GenerateConfig) -> None: self._scale_fragment_detection: float = 1.25 self._scale_minimal_distance: float = 0.8 self._contract_coords: bool = False - self._fixed_charge: int | None = None + self._set_molecular_charge: bool = False + self._molecular_charge: int = 0 def get_identifier(self) -> str: return "generate" @@ -404,20 +405,36 @@ def contract_coords(self, contract_coords: bool): self._contract_coords = contract_coords @property - def fixed_charge(self): + def set_molecular_charge(self): """ - Get the fixed_charge. + Get the set_molecular_charge flag. """ - return self._fixed_charge + return self._set_molecular_charge - @fixed_charge.setter - def fixed_charge(self, fixed_charge: int): + @set_molecular_charge.setter + def set_molecular_charge(self, set_molecular_charge: bool): """ - Set the fixed_charge. + Set the set_molecular_charge flag. """ - if not isinstance(fixed_charge, int): + if not isinstance(set_molecular_charge, bool): + raise TypeError("Set molecular charge should be a boolean.") + self._set_molecular_charge = set_molecular_charge + + @property + def molecular_charge(self): + """ + Get the molecular_charge. + """ + return self._molecular_charge + + @molecular_charge.setter + def molecular_charge(self, molecular_charge: int): + """ + Set the molecular_charge. + """ + if not isinstance(molecular_charge, int): raise TypeError("Fixed charge should be an integer.") - self._fixed_charge = fixed_charge + self._molecular_charge = molecular_charge class RefineConfig(BaseConfig): diff --git a/src/mindlessgen/qm/xtb.py b/src/mindlessgen/qm/xtb.py index fd4411d..e7c1c14 100644 --- a/src/mindlessgen/qm/xtb.py +++ b/src/mindlessgen/qm/xtb.py @@ -241,7 +241,7 @@ def get_xtb_path(binary_name: str | Path | None = None) -> Path: 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. + Check if the remaning number of ligand electrons is even. """ nel = 0 f_electrons = 0 @@ -262,7 +262,6 @@ def check_ligand_uhf(ati: np.ndarray, charge: int) -> None: f_electrons += 102 - atom ln_protons += atom - 3 + 1 ligand_protons = nel - ln_protons - charge - print("ligand protons are:", ligand_protons) # Check if the number of the remaning electrons is even if not ligand_protons % 2 == 0: diff --git a/test/test_generate/test_generate_molecule.py b/test/test_generate/test_generate_molecule.py index 1c8e586..785e219 100644 --- a/test/test_generate/test_generate_molecule.py +++ b/test/test_generate/test_generate_molecule.py @@ -411,8 +411,9 @@ def test_hydrogen_addition( def test_fixed_charge( default_generate_config, ): - """Test the right assinged charge when a fixed charge is given""" - default_generate_config.fixed_charge = 3 + """Test the right assinged charge when a molecular charge is given""" + default_generate_config.molecular_charge = 3 + default_generate_config.set_molecular_charge = True default_generate_config.min_num_atoms = 5 default_generate_config.max_num_atoms = 15 @@ -425,7 +426,8 @@ def test_fixed_charge_and_no_possible_correction( default_generate_config, ): """Test the hydrogen correction for a fixed charge""" - default_generate_config.fixed_charge = 3 + default_generate_config.molecular_charge = 3 + default_generate_config.set_molecular_charge = True default_generate_config.min_num_atoms = 5 default_generate_config.max_num_atoms = 5 default_generate_config.forbidden_elements = "1-10, 12-*" @@ -434,18 +436,64 @@ def test_fixed_charge_and_no_possible_correction( mol = generate_random_molecule(default_generate_config, verbosity=1) assert mol.charge == 3 assert mol.uhf == 0 + assert mol.num_atoms == 5 assert mol.atlist[0] == 0 assert mol.atlist[10] == 5 def test_fixed_charge_hydrogen_correction(default_generate_config): """Test the hydrogen correction for a fixed charge""" - default_generate_config.fixed_charge = 3 - default_generate_config.min_num_atoms = 5 + default_generate_config.molecular_charge = 3 + default_generate_config.set_molecular_charge = True + default_generate_config.min_num_atoms = 7 + default_generate_config.max_num_atoms = 15 + default_generate_config.element_composition = "B:1-1, Ne:2-2, P:1-1, Cl:1-1" + default_generate_config.forbidden_elements = "2-*" + + # Ensure the charge is correctly set + atom_list = generate_atom_list(default_generate_config, verbosity=1) + assert atom_list[0] % 2 == 0 + + +def test_fixed_charge_no_hydrogen_correction(default_generate_config): + """Test the hydrogen correction for a fixed charge""" + default_generate_config.molecular_charge = 2 + default_generate_config.set_molecular_charge = True + default_generate_config.min_num_atoms = 11 + default_generate_config.max_num_atoms = 15 + default_generate_config.element_composition = "H:4-4, B:1-1, Ne:2-2, P:1-1, Cl:1-1" + default_generate_config.forbidden_elements = "1-10, 12-*" + + # Ensure the charge is correctly set + atom_list = generate_atom_list(default_generate_config, verbosity=1) + assert atom_list[0] == 4 + assert atom_list[10] % 2 != 0 + + +def test_fixed_charge_with_lanthanides(default_generate_config): + """Test the hydrogen correction for a fixed charge""" + default_generate_config.molecular_charge = 3 + default_generate_config.set_molecular_charge = True + default_generate_config.min_num_atoms = 7 + default_generate_config.max_num_atoms = 15 + default_generate_config.element_composition = "B:1-1, Ne:2-2, P:1-1, Cl:1-1, Lr:1-1" + default_generate_config.forbidden_elements = "2-*" + + # Ensure the charge is correctly set + atom_list = generate_atom_list(default_generate_config, verbosity=1) + print(atom_list) + assert atom_list[0] % 2 != 0 + + +def test_fixed_charge_with_lanthanides_2(default_generate_config): + """Test the hydrogen correction for a fixed charge""" + default_generate_config.molecular_charge = 2 + default_generate_config.set_molecular_charge = True + default_generate_config.min_num_atoms = 7 default_generate_config.max_num_atoms = 15 - default_generate_config.fixed_elements = "5:1, 10:2, 15:1, 17:1" - default_generate_config.forbidden_elements = "2-10, 11, 12-*" + default_generate_config.element_composition = "B:1-1, Ne:2-2, P:1-1, Cl:1-1, Lu:1-1" + default_generate_config.forbidden_elements = "1-10, 12-*" # Ensure the charge is correctly set atom_list = generate_atom_list(default_generate_config, verbosity=1) - assert np.sum(atom_list[0]) % 2 != 0 + assert atom_list[10] % 2 == 0 From 3627e3b8245be8a5de2e2b9f6cd5cbbfacb7625e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonathan=20Sch=C3=B6ps?= Date: Fri, 8 Nov 2024 17:14:44 +0100 Subject: [PATCH 06/13] Removed an unnessesary import MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Jonathan Schöps --- test/test_molecules/test_miscellaneous.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/test/test_molecules/test_miscellaneous.py b/test/test_molecules/test_miscellaneous.py index 190c51e..b31ce56 100644 --- a/test/test_molecules/test_miscellaneous.py +++ b/test/test_molecules/test_miscellaneous.py @@ -4,11 +4,8 @@ import pytest import numpy as np -from mindlessgen.prog import GenerateConfig # type: ignore from mindlessgen.molecules.miscellaneous import set_random_charge # type: ignore -cfg = GenerateConfig() - # CAUTION: We use 0-based indexing for atoms and molecules! @pytest.mark.parametrize( From d67671eea9bb0685a917dc6c63c3d034ba7d3c09 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonathan=20Sch=C3=B6ps?= Date: Mon, 11 Nov 2024 09:37:26 +0100 Subject: [PATCH 07/13] more convenient function name MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Jonathan Schöps --- src/mindlessgen/molecules/generate_molecule.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/mindlessgen/molecules/generate_molecule.py b/src/mindlessgen/molecules/generate_molecule.py index c9fa746..5bb1d62 100644 --- a/src/mindlessgen/molecules/generate_molecule.py +++ b/src/mindlessgen/molecules/generate_molecule.py @@ -344,11 +344,11 @@ def check_composition(): check_composition() # Check if the number of atoms is within the defined limits check_min_max_atoms() - # If f block elements are included, correct the number of electrons + # If the molecule is not closed shell, add an atom to ensure a closed shell system if cfg.set_molecular_charge: protons, num_atoms = calculate_protons(natoms) nel = protons - cfg.molecular_charge - natoms = fixed_charge_f_block_correction( + natoms = fixed_charge_correction( cfg, natoms, num_atoms, nel, valid_elems, verbosity ) ### ACTUAL WORKFLOW END ### @@ -458,7 +458,7 @@ def check_distances( return True -def fixed_charge_f_block_correction( +def fixed_charge_correction( cfg: GenerateConfig, natoms: np.ndarray, num_atoms: int, @@ -467,7 +467,7 @@ def fixed_charge_f_block_correction( verbosity: int, ) -> np.ndarray: """ - Correct the number of electrons if f block elements are included. + Correct the number of electrons if a fixed charge is given and the molecule is not closed shell. """ f_elem, ligand_protons = calculate_ligand_protons(natoms, nel) # If f block elements are included, cerrect only if the remaning ligand protons are uneven From e7679d70aaf67f462f5290234c12a5cd5179c99e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonathan=20Sch=C3=B6ps?= Date: Thu, 14 Nov 2024 16:44:58 +0100 Subject: [PATCH 08/13] Worked in the suggested changes MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Jonathan Schöps --- mindlessgen.toml | 5 +- src/mindlessgen/cli/cli_parser.py | 7 -- .../molecules/generate_molecule.py | 110 ++++++------------ src/mindlessgen/molecules/miscellaneous.py | 40 ++----- src/mindlessgen/molecules/refinement.py | 29 +++-- src/mindlessgen/prog/config.py | 27 ++--- src/mindlessgen/qm/xtb.py | 59 +++++----- test/test_generate/test_generate_molecule.py | 41 +++++-- test/test_molecules/test_refinement.py | 2 +- 9 files changed, 134 insertions(+), 186 deletions(-) diff --git a/mindlessgen.toml b/mindlessgen.toml index d95025e..6e84c69 100644 --- a/mindlessgen.toml +++ b/mindlessgen.toml @@ -42,9 +42,8 @@ element_composition = "C:2-3, H:1-2, O:1-2, N:1-*" # > CAUTION: This option is overridden by the 'element_composition' option. # > I.e., if an element is specified in 'element_composition' with an occurrence > 0, it will be added to the molecule anyway. forbidden_elements = "57-71, 81-*" -# > Set a charge for the molecule generation. Options: , -set_molecular_charge = false -molecular_charge = 0 +# > Set a charge for the molecule generation. Options: ""(turn it off), "int" +molecular_charge = "" [refine] # > Maximum number of fragment optimization cycles. Options: diff --git a/src/mindlessgen/cli/cli_parser.py b/src/mindlessgen/cli/cli_parser.py index f2517db..65e7ade 100644 --- a/src/mindlessgen/cli/cli_parser.py +++ b/src/mindlessgen/cli/cli_parser.py @@ -143,12 +143,6 @@ def cli_parser(argv: Sequence[str] | None = None) -> dict: required=False, help="Contract the coordinates of the molecule after the coordinats generation.", ) - parser.add_argument( - "--set-charge", - type=bool, - required=False, - help="Set the charge of the molecule.", - ) parser.add_argument( "--molecular-charge", type=int, @@ -292,7 +286,6 @@ def cli_parser(argv: Sequence[str] | None = None) -> dict: "scale_fragment_detection": args_dict["scale_fragment_detection"], "scale_minimal_distance": args_dict["scale_minimal_distance"], "contract_coords": args_dict["contract_coords"], - "set_molecular_charge": args_dict["set_charge"], "molecular_charge": args_dict["molecular_charge"], } # XTB specific arguments diff --git a/src/mindlessgen/molecules/generate_molecule.py b/src/mindlessgen/molecules/generate_molecule.py index 5bb1d62..20d4216 100644 --- a/src/mindlessgen/molecules/generate_molecule.py +++ b/src/mindlessgen/molecules/generate_molecule.py @@ -16,6 +16,7 @@ get_four_d_metals, get_five_d_metals, get_lanthanides, + get_actinides, calculate_ligand_protons, ) @@ -49,7 +50,7 @@ def generate_random_molecule( ati=mol.ati, scale_minimal_distance=config_generate.scale_minimal_distance, ) - if config_generate.set_molecular_charge: + if config_generate.molecular_charge is not None: mol.charge = config_generate.molecular_charge mol.uhf = 0 else: @@ -139,17 +140,16 @@ def generate_atom_list(cfg: GenerateConfig, verbosity: int = 1) -> np.ndarray: + f"Returning: \n{natoms}\n" ) # If the molecular charge is defined, and a fixed element composition is defined, check if the electrons are even. If not raise an error. - if cfg.set_molecular_charge: - protons, num_atoms = calculate_protons(natoms) + if cfg.molecular_charge: + protons = calculate_protons(natoms) nel = protons - cfg.molecular_charge - f_elem, ligand_protons = calculate_ligand_protons(natoms, nel) - if f_elem and ligand_protons % 2 != 0 and fixed_elem: - raise SystemExit( - "Both fixed charge and fixed composition are defined. " - + "Please only define one of them." - + "Or ensure that the fixed composition is closed shell." - ) - elif nel % 2 != 0 and fixed_elem: + f_elem = any( + count > 0 and (i in get_lanthanides() or i in get_actinides()) + for i, count in enumerate(natoms) + ) + if (f_elem and calculate_ligand_protons(natoms, nel) % 2 != 0) or ( + not f_elem and nel % 2 != 0 + ): raise SystemExit( "Both fixed charge and fixed composition are defined. " + "Please only define one of them." @@ -345,12 +345,10 @@ def check_composition(): # Check if the number of atoms is within the defined limits check_min_max_atoms() # If the molecule is not closed shell, add an atom to ensure a closed shell system - if cfg.set_molecular_charge: - protons, num_atoms = calculate_protons(natoms) + if cfg.molecular_charge is not None: + protons = calculate_protons(natoms) nel = protons - cfg.molecular_charge - natoms = fixed_charge_correction( - cfg, natoms, num_atoms, nel, valid_elems, verbosity - ) + natoms = fixed_charge_correction(cfg, natoms, nel, valid_elems, verbosity) ### ACTUAL WORKFLOW END ### return natoms @@ -461,7 +459,6 @@ def check_distances( def fixed_charge_correction( cfg: GenerateConfig, natoms: np.ndarray, - num_atoms: int, nel: int, valid_elems: set[int], verbosity: int, @@ -469,81 +466,45 @@ def fixed_charge_correction( """ Correct the number of electrons if a fixed charge is given and the molecule is not closed shell. """ - f_elem, ligand_protons = calculate_ligand_protons(natoms, nel) - # If f block elements are included, cerrect only if the remaning ligand protons are uneven - if f_elem and ligand_protons % 2 != 0: - if 0 in valid_elems: - natoms = fixed_charge_hydrogen_correction( - cfg, natoms, num_atoms, valid_elems, verbosity - ) - return natoms - else: - natoms = fixed_charge_elem_correction( - cfg, natoms, num_atoms, valid_elems, verbosity - ) + f_elem = any( + count > 0 and (i in get_lanthanides() or i in get_actinides()) + for i, count in enumerate(natoms) + ) + if f_elem: + ligand_protons = calculate_ligand_protons(natoms, nel) + # If f block elements are included, correct only if the remaning ligand protons are uneven + if ligand_protons % 2 != 0: + natoms = fixed_charge_elem_correction(cfg, natoms, valid_elems, verbosity) return natoms # If f block elements are not included, correct if the number of electrons is uneven elif nel % 2 != 0 and f_elem is False: - if 0 in valid_elems: - natoms = fixed_charge_hydrogen_correction( - cfg, natoms, num_atoms, valid_elems, verbosity - ) - return natoms - else: - natoms = fixed_charge_elem_correction( - cfg, natoms, num_atoms, valid_elems, verbosity - ) - return natoms - return natoms - - -def fixed_charge_hydrogen_correction( - cfg: GenerateConfig, - natoms: np.ndarray, - num_atoms: int, - valid_elems: set[int], - verbosity: int, -) -> np.ndarray: - """ - Correct the number of hydrogen atoms if the number of electrons is odd. - """ - min_count, max_count = cfg.element_composition.get(0, (0, cfg.max_num_atoms)) - if max_count is None: - max_count = cfg.max_num_atoms - # Check if adding or removing hydrogen atoms is possible - if natoms[0] < max_count and num_atoms < cfg.max_num_atoms: - natoms[0] += 1 - if verbosity > 1: - print("Adding hydrogen atom...") - return natoms - elif natoms[0] > min_count and num_atoms > cfg.min_num_atoms: - natoms[0] -= 1 - if verbosity > 1: - print("Removing hydrogen atom...") - return natoms - else: - natoms = fixed_charge_elem_correction( - cfg, natoms, num_atoms, valid_elems, verbosity - ) + natoms = fixed_charge_elem_correction(cfg, natoms, valid_elems, verbosity) return natoms + return natoms def fixed_charge_elem_correction( cfg: GenerateConfig, natoms: np.ndarray, - num_atoms: int, valid_elems: set[int], verbosity: int, ) -> np.ndarray: """ Correct the number of atoms if the number of electrons is odd and hydrogen can not be added. """ + num_atoms = np.sum(natoms) # All other elements - odd_atoms = [elem for elem in valid_elems if elem % 2 == 0 and elem != 0] - for random_elem in np.random.permutation(odd_atoms): + random_odd_atoms = np.random.permutation( + [elem for elem in valid_elems if elem % 2 == 0 and elem != 0] + ) + if 0 in valid_elems: + random_odd_atoms = np.insert(random_odd_atoms, 0, 0) + for random_elem in random_odd_atoms: min_count, max_count = cfg.element_composition.get( random_elem, (0, cfg.max_num_atoms) ) + if min_count is None: + min_count = 0 if max_count is None: max_count = cfg.max_num_atoms # Check if adding or removing the random element is possible @@ -557,5 +518,4 @@ def fixed_charge_elem_correction( if verbosity > 1: print(f"Removing atom type {random_elem} for charge...") return natoms - raise RuntimeError("Could not correct the odd number of electrons.") - return natoms + raise RuntimeError("Could not correct the odd number of electrons.") diff --git a/src/mindlessgen/molecules/miscellaneous.py b/src/mindlessgen/molecules/miscellaneous.py index 8d3747e..d205ba9 100644 --- a/src/mindlessgen/molecules/miscellaneous.py +++ b/src/mindlessgen/molecules/miscellaneous.py @@ -55,7 +55,7 @@ def set_random_charge( f"Number of protons from ligands (assuming negative charge): {ligand_protons}" ) - elif ligand_protons % 2 == 0: + if ligand_protons % 2 == 0: charge = 0 else: charge = 1 @@ -84,43 +84,27 @@ def set_random_charge( return charge, uhf -def calculate_protons(natoms: np.ndarray) -> tuple[int, int]: +def calculate_protons(natoms: np.ndarray) -> int: """ Calculate the number of protons in a molecule from the atom list. """ protons = 0 - num_atoms = 0 - temp_count = 0 - for atom in natoms: - temp_count += 1 - protons += atom * temp_count - if atom > 0: - num_atoms += atom - return protons, num_atoms + for i, atom_count in enumerate(natoms): + protons += atom_count * (i + 1) + return protons -def calculate_ligand_protons(natoms: np.ndarray, nel) -> tuple[bool, int]: +def calculate_ligand_protons(natoms: np.ndarray, nel: int) -> int: """ Calculate the number of ligand protons in a molecule if lanthanides or actinides are within the molecule. """ - ligand_protons = 0 - f_protons = 0 - f_elem = False - temp_count = -1 - for atom in natoms: - temp_count += 1 - if ( - atom > 0 - and temp_count in get_lanthanides() - or atom > 0 - and temp_count in get_actinides() - ): - f_protons += atom * ( - temp_count - 3 + 1 - ) # subtract 3 to get the number of protons in the Ln3+ ion - f_elem = True + f_protons = sum( + atom * (i - 2) + for i, atom in enumerate(natoms) + if (i in get_lanthanides() or i in get_actinides()) + ) ligand_protons = nel - f_protons - return f_elem, ligand_protons + return ligand_protons def get_alkali_metals() -> list[int]: diff --git a/src/mindlessgen/molecules/refinement.py b/src/mindlessgen/molecules/refinement.py index fc829c3..036f813 100644 --- a/src/mindlessgen/molecules/refinement.py +++ b/src/mindlessgen/molecules/refinement.py @@ -13,6 +13,8 @@ set_random_charge, calculate_protons, calculate_ligand_protons, + get_lanthanides, + get_actinides, ) COV_RADII = "pyykko" @@ -54,7 +56,7 @@ def iterative_optimization( # Detect fragments from the optimized molecule fragmols = detect_fragments( mol=rev_mol, - config_generate=config_generate, + molecular_charge=config_generate.molecular_charge, vdw_scaling=config_generate.scale_fragment_detection, verbosity=verbosity, ) @@ -114,14 +116,19 @@ def iterative_optimization( f"Element {ati} is overrepresented " + f"in the largest fragment in cycle {cycle + 1}." ) - if config_generate.set_molecular_charge: - protons, num_atoms = calculate_protons(fragmols[0].atlist) + if config_generate.molecular_charge is not None: + protons = calculate_protons(fragmols[0].atlist) nel = protons - config_generate.molecular_charge - f_elem, ligand_protons = calculate_ligand_protons(fragmols[0].atlist, nel) - if ligand_protons % 2 != 0: - raise RuntimeError( - f"Number of electrons in the largest fragment in cycle {cycle + 1} is odd." - ) + f_elem = any( + count > 0 and (i in get_lanthanides() or i in get_actinides()) + for i, count in enumerate(fragmols[0].atlist) + ) + if f_elem: + ligand_protons = calculate_ligand_protons(fragmols[0].atlist, nel) + if ligand_protons % 2 != 0: + raise RuntimeError( + f"Number of electrons in the largest fragment in cycle {cycle + 1} is odd." + ) elif nel % 2 != 0: raise RuntimeError( f"Number of electrons in the largest fragment in cycle {cycle + 1} is odd." @@ -153,7 +160,7 @@ def iterative_optimization( def detect_fragments( mol: Molecule, - config_generate: GenerateConfig, + molecular_charge: int, vdw_scaling: float, verbosity: int = 1, ) -> list[Molecule]: @@ -220,8 +227,8 @@ def detect_fragments( for atom in fragment_molecule.ati: fragment_molecule.atlist[atom] += 1 # Update the charge of the fragment molecule - if config_generate.set_molecular_charge: - fragment_molecule.charge = config_generate.molecular_charge + if molecular_charge is not None: + fragment_molecule.charge = molecular_charge fragment_molecule.uhf = 0 else: fragment_molecule.charge, fragment_molecule.uhf = set_random_charge( diff --git a/src/mindlessgen/prog/config.py b/src/mindlessgen/prog/config.py index 6ec54e6..3bf3794 100644 --- a/src/mindlessgen/prog/config.py +++ b/src/mindlessgen/prog/config.py @@ -183,8 +183,7 @@ def __init__(self: GenerateConfig) -> None: self._scale_fragment_detection: float = 1.25 self._scale_minimal_distance: float = 0.8 self._contract_coords: bool = True - self._set_molecular_charge: bool = False - self._molecular_charge: int = 0 + self._molecular_charge: int | None = None def get_identifier(self) -> str: return "generate" @@ -404,22 +403,6 @@ def contract_coords(self, contract_coords: bool): raise TypeError("Contract coords should be a boolean.") self._contract_coords = contract_coords - @property - def set_molecular_charge(self): - """ - Get the set_molecular_charge flag. - """ - return self._set_molecular_charge - - @set_molecular_charge.setter - def set_molecular_charge(self, set_molecular_charge: bool): - """ - Set the set_molecular_charge flag. - """ - if not isinstance(set_molecular_charge, bool): - raise TypeError("Set molecular charge should be a boolean.") - self._set_molecular_charge = set_molecular_charge - @property def molecular_charge(self): """ @@ -428,11 +411,15 @@ def molecular_charge(self): return self._molecular_charge @molecular_charge.setter - def molecular_charge(self, molecular_charge: int): + def molecular_charge(self, molecular_charge_str: str): """ Set the molecular_charge. """ - if not isinstance(molecular_charge, int): + if molecular_charge_str.strip() == "": + return None + else: + molecular_charge = int(molecular_charge_str) + if not isinstance(molecular_charge, int | None): raise TypeError("Fixed charge should be an integer.") self._molecular_charge = molecular_charge diff --git a/src/mindlessgen/qm/xtb.py b/src/mindlessgen/qm/xtb.py index e7c1c14..efa563a 100644 --- a/src/mindlessgen/qm/xtb.py +++ b/src/mindlessgen/qm/xtb.py @@ -10,7 +10,12 @@ from ..molecules import Molecule from .base import QMMethod from ..prog import XTBConfig -from ..molecules.miscellaneous import get_lanthanides, get_actinides +from ..molecules.miscellaneous import ( + calculate_ligand_protons, + calculate_protons, + get_lanthanides, + get_actinides, +) class XTB(QMMethod): @@ -41,8 +46,11 @@ def optimize( if np.any(molecule.ati > 85): super_heavy_elements = True molecule.ati[molecule.ati > 85] -= 32 - if np.any(np.isin(molecule.ati, get_lanthanides())): - check_ligand_uhf(molecule.ati, molecule.charge) + if np.any( + np.isin(molecule.ati, get_lanthanides()) + or np.isin(molecule.ati, get_actinides()) + ): + check_ligand_uhf(molecule.atlist, 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. @@ -84,7 +92,10 @@ 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())): + if np.any( + np.isin(molecule.ati, get_lanthanides()) + or np.isin(molecule.ati, get_actinides()) + ): # Reset the UHF value to the original value before returning the optimized molecule. optimized_molecule.uhf = uhf_original if super_heavy_elements: @@ -102,8 +113,11 @@ def singlepoint(self, molecule: Molecule, verbosity: int = 1) -> str: if np.any(molecule.ati > 85): super_heavy_elements = True molecule.ati[molecule.ati > 85] -= 32 - if np.any(np.isin(molecule.ati, get_lanthanides())): - check_ligand_uhf(molecule.ati, molecule.charge) + if np.any( + np.isin(molecule.ati, get_lanthanides()) + or np.isin(molecule.ati, get_actinides()) + ): + check_ligand_uhf(molecule.atlist, 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. @@ -140,7 +154,10 @@ 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())): + if np.any( + np.isin(molecule.ati, get_lanthanides()) + or np.isin(molecule.ati, get_actinides()) + ): molecule.uhf = uhf_original if super_heavy_elements: # Reset the atomic numbers to the original values before returning the optimized molecule. @@ -239,32 +256,14 @@ def get_xtb_path(binary_name: str | Path | None = None) -> Path: raise ImportError("'xtb' binary could not be found.") -def check_ligand_uhf(ati: np.ndarray, charge: int) -> None: +def check_ligand_uhf(atlist: np.ndarray, charge: int) -> None: """ Check if the remaning number of ligand electrons is even. """ - nel = 0 - f_electrons = 0 - ligand_protons = 0 - ln_protons = 0 - for atom in ati: - nel += atom + 1 - if atom in get_lanthanides(): - if atom < 64: - f_electrons += atom - 56 - else: - f_electrons += 70 - atom - ln_protons += atom - 3 + 1 - if atom in get_actinides(): - if atom < 96: - f_electrons += atom - 88 - else: - f_electrons += 102 - atom - ln_protons += atom - 3 + 1 - ligand_protons = nel - ln_protons - charge - - # Check if the number of the remaning electrons is even - if not ligand_protons % 2 == 0: + protons = calculate_protons(atlist) + nel = protons - charge + ligand_protons = calculate_ligand_protons(atlist, nel) + if ligand_protons % 2 != 0: raise SystemExit( "The number of electrons in the ligand is not even. Please check the input." ) diff --git a/test/test_generate/test_generate_molecule.py b/test/test_generate/test_generate_molecule.py index 785e219..7369e1a 100644 --- a/test/test_generate/test_generate_molecule.py +++ b/test/test_generate/test_generate_molecule.py @@ -412,7 +412,7 @@ def test_fixed_charge( default_generate_config, ): """Test the right assinged charge when a molecular charge is given""" - default_generate_config.molecular_charge = 3 + default_generate_config.molecular_charge = "3" default_generate_config.set_molecular_charge = True default_generate_config.min_num_atoms = 5 default_generate_config.max_num_atoms = 15 @@ -420,13 +420,14 @@ def test_fixed_charge( # Ensure the charge is correctly set mol = generate_random_molecule(default_generate_config, verbosity=1) assert mol.charge == 3 + assert mol.uhf == 0 def test_fixed_charge_and_no_possible_correction( default_generate_config, ): """Test the hydrogen correction for a fixed charge""" - default_generate_config.molecular_charge = 3 + default_generate_config.molecular_charge = "3" default_generate_config.set_molecular_charge = True default_generate_config.min_num_atoms = 5 default_generate_config.max_num_atoms = 5 @@ -441,30 +442,48 @@ def test_fixed_charge_and_no_possible_correction( assert mol.atlist[10] == 5 +def test_fixed_charge_and_fixed_composition( + default_generate_config, +): + """Test the hydrogen correction for a fixed charge""" + default_generate_config.molecular_charge = "3" + default_generate_config.set_molecular_charge = True + default_generate_config.min_num_atoms = 5 + default_generate_config.max_num_atoms = 10 + default_generate_config.element_composition = "H:5-5, C:2-2, N:1-1, O:1-1" + + # Check if the right system exit is raised + with pytest.raises( + SystemExit, + match="Both fixed charge and fixed composition are defined. Please only define one of them.Or ensure that the fixed composition is closed shell.", + ): + generate_random_molecule(default_generate_config, verbosity=1) + + def test_fixed_charge_hydrogen_correction(default_generate_config): """Test the hydrogen correction for a fixed charge""" - default_generate_config.molecular_charge = 3 + default_generate_config.molecular_charge = "3" default_generate_config.set_molecular_charge = True default_generate_config.min_num_atoms = 7 default_generate_config.max_num_atoms = 15 default_generate_config.element_composition = "B:1-1, Ne:2-2, P:1-1, Cl:1-1" default_generate_config.forbidden_elements = "2-*" - # Ensure the charge is correctly set + # Ensure the right hydrogen correction is applied atom_list = generate_atom_list(default_generate_config, verbosity=1) assert atom_list[0] % 2 == 0 def test_fixed_charge_no_hydrogen_correction(default_generate_config): """Test the hydrogen correction for a fixed charge""" - default_generate_config.molecular_charge = 2 + default_generate_config.molecular_charge = "2" default_generate_config.set_molecular_charge = True default_generate_config.min_num_atoms = 11 default_generate_config.max_num_atoms = 15 default_generate_config.element_composition = "H:4-4, B:1-1, Ne:2-2, P:1-1, Cl:1-1" default_generate_config.forbidden_elements = "1-10, 12-*" - # Ensure the charge is correctly set + # Ensure the right atom correction is applied atom_list = generate_atom_list(default_generate_config, verbosity=1) assert atom_list[0] == 4 assert atom_list[10] % 2 != 0 @@ -472,28 +491,28 @@ def test_fixed_charge_no_hydrogen_correction(default_generate_config): def test_fixed_charge_with_lanthanides(default_generate_config): """Test the hydrogen correction for a fixed charge""" - default_generate_config.molecular_charge = 3 + default_generate_config.molecular_charge = "3" default_generate_config.set_molecular_charge = True default_generate_config.min_num_atoms = 7 default_generate_config.max_num_atoms = 15 default_generate_config.element_composition = "B:1-1, Ne:2-2, P:1-1, Cl:1-1, Lr:1-1" default_generate_config.forbidden_elements = "2-*" - # Ensure the charge is correctly set + # Ensure the right ligand correction is applied atom_list = generate_atom_list(default_generate_config, verbosity=1) - print(atom_list) assert atom_list[0] % 2 != 0 def test_fixed_charge_with_lanthanides_2(default_generate_config): """Test the hydrogen correction for a fixed charge""" - default_generate_config.molecular_charge = 2 + default_generate_config.molecular_charge = "2" default_generate_config.set_molecular_charge = True default_generate_config.min_num_atoms = 7 default_generate_config.max_num_atoms = 15 default_generate_config.element_composition = "B:1-1, Ne:2-2, P:1-1, Cl:1-1, Lu:1-1" default_generate_config.forbidden_elements = "1-10, 12-*" - # Ensure the charge is correctly set + # Ensure the right ligand correction is applied atom_list = generate_atom_list(default_generate_config, verbosity=1) + assert atom_list[0] == 0 assert atom_list[10] % 2 == 0 diff --git a/test/test_molecules/test_refinement.py b/test/test_molecules/test_refinement.py index 087105f..4011e3c 100644 --- a/test/test_molecules/test_refinement.py +++ b/test/test_molecules/test_refinement.py @@ -108,7 +108,7 @@ def test_detect_fragments_H6O2B2Ne2I1Os1Tl1( """ fragmols = detect_fragments( mol=mol_H6O2B2Ne2I1Os1Tl1, - config_generate=GenerateConfig(), + molecular_charge=GenerateConfig().molecular_charge, vdw_scaling=1.3333, verbosity=0, ) From 9c6903b0f6b939c8ce73ec45d81ce68df5822365 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonathan=20Sch=C3=B6ps?= Date: Thu, 14 Nov 2024 16:55:08 +0100 Subject: [PATCH 09/13] Added a changelog entry MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Jonathan Schöps --- CHANGELOG.md | 1 + src/mindlessgen/qm/xtb.py | 37 ++++++++----------------------------- 2 files changed, 9 insertions(+), 29 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1cf4174..7537601 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -28,6 +28,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - function which is able to printout the xyz coordinates to the terminal similar to the `.xyz` layout - elements 87 to 103 are accessible via the element composition. If `xtb` is the engine, the elements will be replaced by their lighter homologues. - support for `python-3.13` +- function which is abel to set a molecular charge and ensure `uhf = 0` ### Breaking Changes - Removal of the `dist_threshold` flag and in the `-toml` file. diff --git a/src/mindlessgen/qm/xtb.py b/src/mindlessgen/qm/xtb.py index efa563a..9431e21 100644 --- a/src/mindlessgen/qm/xtb.py +++ b/src/mindlessgen/qm/xtb.py @@ -11,8 +11,6 @@ from .base import QMMethod from ..prog import XTBConfig from ..molecules.miscellaneous import ( - calculate_ligand_protons, - calculate_protons, get_lanthanides, get_actinides, ) @@ -46,11 +44,9 @@ def optimize( if np.any(molecule.ati > 85): super_heavy_elements = True molecule.ati[molecule.ati > 85] -= 32 - if np.any( - np.isin(molecule.ati, get_lanthanides()) - or np.isin(molecule.ati, get_actinides()) + if np.any(np.isin(molecule.ati, get_lanthanides())) or np.any( + np.isin(molecule.ati, get_actinides()) ): - check_ligand_uhf(molecule.atlist, 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. @@ -92,9 +88,8 @@ 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()) - or np.isin(molecule.ati, get_actinides()) + if np.any(np.isin(molecule.ati, get_lanthanides())) or np.any( + np.isin(molecule.ati, get_actinides()) ): # Reset the UHF value to the original value before returning the optimized molecule. optimized_molecule.uhf = uhf_original @@ -113,11 +108,9 @@ def singlepoint(self, molecule: Molecule, verbosity: int = 1) -> str: if np.any(molecule.ati > 85): super_heavy_elements = True molecule.ati[molecule.ati > 85] -= 32 - if np.any( - np.isin(molecule.ati, get_lanthanides()) - or np.isin(molecule.ati, get_actinides()) + if np.any(np.isin(molecule.ati, get_lanthanides())) or np.any( + np.isin(molecule.ati, get_actinides()) ): - check_ligand_uhf(molecule.atlist, 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. @@ -154,9 +147,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()) - or np.isin(molecule.ati, get_actinides()) + if np.any(np.isin(molecule.ati, get_lanthanides())) or np.any( + np.isin(molecule.ati, get_actinides()) ): molecule.uhf = uhf_original if super_heavy_elements: @@ -254,16 +246,3 @@ 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(atlist: np.ndarray, charge: int) -> None: - """ - Check if the remaning number of ligand electrons is even. - """ - protons = calculate_protons(atlist) - nel = protons - charge - ligand_protons = calculate_ligand_protons(atlist, nel) - if ligand_protons % 2 != 0: - raise SystemExit( - "The number of electrons in the ligand is not even. Please check the input." - ) From d8be3d0d9f2d3417dd2f48dca7e14ab5d80c116f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonathan=20Sch=C3=B6ps?= Date: Thu, 14 Nov 2024 17:35:17 +0100 Subject: [PATCH 10/13] Changes a bug with an array type MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Jonathan Schöps --- src/mindlessgen/molecules/generate_molecule.py | 9 ++++++++- test/test_generate/test_generate_molecule.py | 7 ------- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/mindlessgen/molecules/generate_molecule.py b/src/mindlessgen/molecules/generate_molecule.py index 04a6f0a..dbdfc3c 100644 --- a/src/mindlessgen/molecules/generate_molecule.py +++ b/src/mindlessgen/molecules/generate_molecule.py @@ -499,11 +499,17 @@ def fixed_charge_elem_correction( """ num_atoms = np.sum(natoms) # All other elements + random_odd_atoms = np.array([], dtype=int) + print(f"Odd atoms: {random_odd_atoms.dtype}") random_odd_atoms = np.random.permutation( [elem for elem in valid_elems if elem % 2 == 0 and elem != 0] ) + random_odd_atoms = np.array(random_odd_atoms, dtype=int) + print(f"Valid elements: {valid_elems}, random odd atoms: {random_odd_atoms.dtype}") + print(f"Valid elements: {random_odd_atoms}") if 0 in valid_elems: - random_odd_atoms = np.insert(random_odd_atoms, 0, 0) + random_odd_atoms = np.insert(random_odd_atoms, 0, int(0)) + print(f"Random odd atoms: {random_odd_atoms}") for random_elem in random_odd_atoms: min_count, max_count = cfg.element_composition.get( random_elem, (0, cfg.max_num_atoms) @@ -512,6 +518,7 @@ def fixed_charge_elem_correction( min_count = 0 if max_count is None: max_count = cfg.max_num_atoms + print(f"Random element: {random_elem}") # Check if adding or removing the random element is possible if natoms[random_elem] < max_count and num_atoms < cfg.max_num_atoms: natoms[random_elem] += 1 diff --git a/test/test_generate/test_generate_molecule.py b/test/test_generate/test_generate_molecule.py index 7369e1a..c7f0ca1 100644 --- a/test/test_generate/test_generate_molecule.py +++ b/test/test_generate/test_generate_molecule.py @@ -413,7 +413,6 @@ def test_fixed_charge( ): """Test the right assinged charge when a molecular charge is given""" default_generate_config.molecular_charge = "3" - default_generate_config.set_molecular_charge = True default_generate_config.min_num_atoms = 5 default_generate_config.max_num_atoms = 15 @@ -428,7 +427,6 @@ def test_fixed_charge_and_no_possible_correction( ): """Test the hydrogen correction for a fixed charge""" default_generate_config.molecular_charge = "3" - default_generate_config.set_molecular_charge = True default_generate_config.min_num_atoms = 5 default_generate_config.max_num_atoms = 5 default_generate_config.forbidden_elements = "1-10, 12-*" @@ -447,7 +445,6 @@ def test_fixed_charge_and_fixed_composition( ): """Test the hydrogen correction for a fixed charge""" default_generate_config.molecular_charge = "3" - default_generate_config.set_molecular_charge = True default_generate_config.min_num_atoms = 5 default_generate_config.max_num_atoms = 10 default_generate_config.element_composition = "H:5-5, C:2-2, N:1-1, O:1-1" @@ -463,7 +460,6 @@ def test_fixed_charge_and_fixed_composition( def test_fixed_charge_hydrogen_correction(default_generate_config): """Test the hydrogen correction for a fixed charge""" default_generate_config.molecular_charge = "3" - default_generate_config.set_molecular_charge = True default_generate_config.min_num_atoms = 7 default_generate_config.max_num_atoms = 15 default_generate_config.element_composition = "B:1-1, Ne:2-2, P:1-1, Cl:1-1" @@ -477,7 +473,6 @@ def test_fixed_charge_hydrogen_correction(default_generate_config): def test_fixed_charge_no_hydrogen_correction(default_generate_config): """Test the hydrogen correction for a fixed charge""" default_generate_config.molecular_charge = "2" - default_generate_config.set_molecular_charge = True default_generate_config.min_num_atoms = 11 default_generate_config.max_num_atoms = 15 default_generate_config.element_composition = "H:4-4, B:1-1, Ne:2-2, P:1-1, Cl:1-1" @@ -492,7 +487,6 @@ def test_fixed_charge_no_hydrogen_correction(default_generate_config): def test_fixed_charge_with_lanthanides(default_generate_config): """Test the hydrogen correction for a fixed charge""" default_generate_config.molecular_charge = "3" - default_generate_config.set_molecular_charge = True default_generate_config.min_num_atoms = 7 default_generate_config.max_num_atoms = 15 default_generate_config.element_composition = "B:1-1, Ne:2-2, P:1-1, Cl:1-1, Lr:1-1" @@ -506,7 +500,6 @@ def test_fixed_charge_with_lanthanides(default_generate_config): def test_fixed_charge_with_lanthanides_2(default_generate_config): """Test the hydrogen correction for a fixed charge""" default_generate_config.molecular_charge = "2" - default_generate_config.set_molecular_charge = True default_generate_config.min_num_atoms = 7 default_generate_config.max_num_atoms = 15 default_generate_config.element_composition = "B:1-1, Ne:2-2, P:1-1, Cl:1-1, Lu:1-1" From 22685531d4f738c70d2231beefef4ac8bf4d88cd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonathan=20Sch=C3=B6ps?= Date: Thu, 14 Nov 2024 17:41:01 +0100 Subject: [PATCH 11/13] removed left over print statements MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Jonathan Schöps --- src/mindlessgen/molecules/generate_molecule.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/mindlessgen/molecules/generate_molecule.py b/src/mindlessgen/molecules/generate_molecule.py index dbdfc3c..a7f6cee 100644 --- a/src/mindlessgen/molecules/generate_molecule.py +++ b/src/mindlessgen/molecules/generate_molecule.py @@ -500,16 +500,12 @@ def fixed_charge_elem_correction( num_atoms = np.sum(natoms) # All other elements random_odd_atoms = np.array([], dtype=int) - print(f"Odd atoms: {random_odd_atoms.dtype}") random_odd_atoms = np.random.permutation( [elem for elem in valid_elems if elem % 2 == 0 and elem != 0] ) random_odd_atoms = np.array(random_odd_atoms, dtype=int) - print(f"Valid elements: {valid_elems}, random odd atoms: {random_odd_atoms.dtype}") - print(f"Valid elements: {random_odd_atoms}") if 0 in valid_elems: random_odd_atoms = np.insert(random_odd_atoms, 0, int(0)) - print(f"Random odd atoms: {random_odd_atoms}") for random_elem in random_odd_atoms: min_count, max_count = cfg.element_composition.get( random_elem, (0, cfg.max_num_atoms) @@ -518,7 +514,6 @@ def fixed_charge_elem_correction( min_count = 0 if max_count is None: max_count = cfg.max_num_atoms - print(f"Random element: {random_elem}") # Check if adding or removing the random element is possible if natoms[random_elem] < max_count and num_atoms < cfg.max_num_atoms: natoms[random_elem] += 1 From 869846a3ebd7afa31ef7c1b350d0a2f53b4f3f07 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonathan=20Sch=C3=B6ps?= Date: Fri, 15 Nov 2024 14:21:38 +0100 Subject: [PATCH 12/13] Implemented requested changes MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Jonathan Schöps --- CHANGELOG.md | 2 +- mindlessgen.toml | 4 ++-- src/mindlessgen/cli/cli_parser.py | 2 +- .../molecules/generate_molecule.py | 23 ++++++++----------- src/mindlessgen/molecules/miscellaneous.py | 19 ++++++++------- src/mindlessgen/molecules/refinement.py | 6 ++--- src/mindlessgen/prog/config.py | 16 +++++++------ src/mindlessgen/qm/xtb.py | 4 ++-- .../test_config/test_config_set_attributes.py | 18 +++++++++++++++ 9 files changed, 57 insertions(+), 37 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7537601..a597057 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -28,7 +28,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - function which is able to printout the xyz coordinates to the terminal similar to the `.xyz` layout - elements 87 to 103 are accessible via the element composition. If `xtb` is the engine, the elements will be replaced by their lighter homologues. - support for `python-3.13` -- function which is abel to set a molecular charge and ensure `uhf = 0` +- option to set a fixed molecular charge, while ensuring `uhf = 0` ### Breaking Changes - Removal of the `dist_threshold` flag and in the `-toml` file. diff --git a/mindlessgen.toml b/mindlessgen.toml index 6e84c69..538e071 100644 --- a/mindlessgen.toml +++ b/mindlessgen.toml @@ -42,8 +42,8 @@ element_composition = "C:2-3, H:1-2, O:1-2, N:1-*" # > CAUTION: This option is overridden by the 'element_composition' option. # > I.e., if an element is specified in 'element_composition' with an occurrence > 0, it will be added to the molecule anyway. forbidden_elements = "57-71, 81-*" -# > Set a charge for the molecule generation. Options: ""(turn it off), "int" -molecular_charge = "" +# > Set a charge for the molecule generation. Options: "none" or "" (turn it off), "int" or +molecular_charge = "none" [refine] # > Maximum number of fragment optimization cycles. Options: diff --git a/src/mindlessgen/cli/cli_parser.py b/src/mindlessgen/cli/cli_parser.py index 65e7ade..e68ef8d 100644 --- a/src/mindlessgen/cli/cli_parser.py +++ b/src/mindlessgen/cli/cli_parser.py @@ -145,7 +145,7 @@ def cli_parser(argv: Sequence[str] | None = None) -> dict: ) parser.add_argument( "--molecular-charge", - type=int, + type=str, required=False, help="Define the charge of the molecule.", ) diff --git a/src/mindlessgen/molecules/generate_molecule.py b/src/mindlessgen/molecules/generate_molecule.py index a7f6cee..cc6408f 100644 --- a/src/mindlessgen/molecules/generate_molecule.py +++ b/src/mindlessgen/molecules/generate_molecule.py @@ -17,7 +17,7 @@ get_five_d_metals, get_lanthanides, get_actinides, - calculate_ligand_protons, + calculate_ligand_electrons, ) @@ -56,7 +56,6 @@ def generate_random_molecule( else: mol.charge, mol.uhf = set_random_charge(mol.ati, verbosity) mol.set_name_from_formula() - if verbosity > 1: print(mol) @@ -150,7 +149,7 @@ def generate_atom_list(cfg: GenerateConfig, verbosity: int = 1) -> np.ndarray: count > 0 and (i in get_lanthanides() or i in get_actinides()) for i, count in enumerate(natoms) ) - if (f_elem and calculate_ligand_protons(natoms, nel) % 2 != 0) or ( + if (f_elem and calculate_ligand_electrons(natoms, nel) % 2 != 0) or ( not f_elem and nel % 2 != 0 ): raise SystemExit( @@ -476,13 +475,13 @@ def fixed_charge_correction( for i, count in enumerate(natoms) ) if f_elem: - ligand_protons = calculate_ligand_protons(natoms, nel) + ligand_electrons = calculate_ligand_electrons(natoms, nel) # If f block elements are included, correct only if the remaning ligand protons are uneven - if ligand_protons % 2 != 0: + if ligand_electrons % 2 != 0: natoms = fixed_charge_elem_correction(cfg, natoms, valid_elems, verbosity) return natoms # If f block elements are not included, correct if the number of electrons is uneven - elif nel % 2 != 0 and f_elem is False: + elif nel % 2 != 0: natoms = fixed_charge_elem_correction(cfg, natoms, valid_elems, verbosity) return natoms return natoms @@ -495,17 +494,15 @@ def fixed_charge_elem_correction( verbosity: int, ) -> np.ndarray: """ - Correct the number of atoms if the number of electrons is odd and hydrogen can not be added. + Correct the number of atoms if the number of electrons is odd and a molecular charge is set. """ num_atoms = np.sum(natoms) # All other elements - random_odd_atoms = np.array([], dtype=int) - random_odd_atoms = np.random.permutation( - [elem for elem in valid_elems if elem % 2 == 0 and elem != 0] - ) - random_odd_atoms = np.array(random_odd_atoms, dtype=int) + rng = np.random.default_rng() + odd_atoms = np.array([elem for elem in valid_elems if elem % 2 == 0], dtype=int) + random_odd_atoms = rng.permutation(odd_atoms) if 0 in valid_elems: - random_odd_atoms = np.insert(random_odd_atoms, 0, int(0)) + random_odd_atoms = np.insert(random_odd_atoms, 0, 0) for random_elem in random_odd_atoms: min_count, max_count = cfg.element_composition.get( random_elem, (0, cfg.max_num_atoms) diff --git a/src/mindlessgen/molecules/miscellaneous.py b/src/mindlessgen/molecules/miscellaneous.py index d4abeff..d0ad736 100644 --- a/src/mindlessgen/molecules/miscellaneous.py +++ b/src/mindlessgen/molecules/miscellaneous.py @@ -95,17 +95,20 @@ def calculate_protons(natoms: np.ndarray) -> int: return protons -def calculate_ligand_protons(natoms: np.ndarray, nel: int) -> int: +def calculate_ligand_electrons(natoms: np.ndarray, nel: int) -> int: """ - Calculate the number of ligand protons in a molecule if lanthanides or actinides are within the molecule. + Calculate the number of ligand electrons in a molecule if lanthanides or actinides are within the molecule. """ - f_protons = sum( - atom * (i - 2) - for i, atom in enumerate(natoms) - if (i in get_lanthanides() or i in get_actinides()) + f_electrons = sum( + occurence + * ( + ati - 3 + 1 + ) # subtract 3 to get the number of electrons for an Ln3+ (Ac3+) ion. + for ati, occurence in enumerate(natoms) + if (ati in get_lanthanides() or ati in get_actinides()) ) - ligand_protons = nel - f_protons - return ligand_protons + ligand_electrons = nel - f_electrons + return ligand_electrons def get_alkali_metals() -> list[int]: diff --git a/src/mindlessgen/molecules/refinement.py b/src/mindlessgen/molecules/refinement.py index 036f813..4e884c9 100644 --- a/src/mindlessgen/molecules/refinement.py +++ b/src/mindlessgen/molecules/refinement.py @@ -12,7 +12,7 @@ from .miscellaneous import ( set_random_charge, calculate_protons, - calculate_ligand_protons, + calculate_ligand_electrons, get_lanthanides, get_actinides, ) @@ -124,8 +124,8 @@ def iterative_optimization( for i, count in enumerate(fragmols[0].atlist) ) if f_elem: - ligand_protons = calculate_ligand_protons(fragmols[0].atlist, nel) - if ligand_protons % 2 != 0: + ligand_electrons = calculate_ligand_electrons(fragmols[0].atlist, nel) + if ligand_electrons % 2 != 0: raise RuntimeError( f"Number of electrons in the largest fragment in cycle {cycle + 1} is odd." ) diff --git a/src/mindlessgen/prog/config.py b/src/mindlessgen/prog/config.py index 3bf3794..1d6aefa 100644 --- a/src/mindlessgen/prog/config.py +++ b/src/mindlessgen/prog/config.py @@ -411,17 +411,19 @@ def molecular_charge(self): return self._molecular_charge @molecular_charge.setter - def molecular_charge(self, molecular_charge_str: str): + def molecular_charge(self, molecular_charge: str | int): """ Set the molecular_charge. """ - if molecular_charge_str.strip() == "": - return None + if isinstance(molecular_charge, str): + if molecular_charge.lower() == "none" or molecular_charge == "": + self._molecular_charge = None + else: + self._molecular_charge = int(molecular_charge) + elif isinstance(molecular_charge, int): + self._molecular_charge = molecular_charge else: - molecular_charge = int(molecular_charge_str) - if not isinstance(molecular_charge, int | None): - raise TypeError("Fixed charge should be an integer.") - self._molecular_charge = molecular_charge + raise TypeError("Molecular charge should be a string or an integer.") class RefineConfig(BaseConfig): diff --git a/src/mindlessgen/qm/xtb.py b/src/mindlessgen/qm/xtb.py index 9431e21..3569d49 100644 --- a/src/mindlessgen/qm/xtb.py +++ b/src/mindlessgen/qm/xtb.py @@ -40,8 +40,8 @@ def optimize( Optimize a molecule using xtb. """ super_heavy_elements = False - ati_original = molecule.ati.copy() if np.any(molecule.ati > 85): + ati_original = molecule.ati.copy() super_heavy_elements = True molecule.ati[molecule.ati > 85] -= 32 if np.any(np.isin(molecule.ati, get_lanthanides())) or np.any( @@ -95,7 +95,7 @@ def optimize( optimized_molecule.uhf = uhf_original if super_heavy_elements: # Reset the atomic numbers to the original values before returning the optimized molecule. - optimized_molecule.ati = ati_original + optimized_molecule.ati = ati_original # pylint: disable=E0606 optimized_molecule.atlist = molecule.atlist return optimized_molecule diff --git a/test/test_config/test_config_set_attributes.py b/test/test_config/test_config_set_attributes.py index feb7ae0..0af85f3 100644 --- a/test/test_config/test_config_set_attributes.py +++ b/test/test_config/test_config_set_attributes.py @@ -123,6 +123,24 @@ def test_generate_config_default_values(property_name, initial_value): assert getattr(config, property_name) == initial_value +# Test for fixed charge +@pytest.mark.parametrize( + "property_name, valid_value, expected_value", + [ + ("molecular_charge", 0, 0), + ("molecular_charge", "1", 1), + ("molecular_charge", "none", None), + ("molecular_charge", "", None), + ], +) +def test_generate_config_molecular_charge(property_name, valid_value, expected_value): + config = GenerateConfig() + + # Test valid value + setattr(config, property_name, valid_value) + assert getattr(config, property_name) == expected_value + + # Tests for RefineConfig @pytest.mark.parametrize( "property_name, valid_value, invalid_value, expected_exception", From ccb96c7aa8fc7f6765fb53baabae49139bb1671d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonathan=20Sch=C3=B6ps?= Date: Fri, 15 Nov 2024 14:58:23 +0100 Subject: [PATCH 13/13] Added a new test to raise a ValueError MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Jonathan Schöps --- mindlessgen.toml | 2 +- src/mindlessgen/molecules/miscellaneous.py | 8 ++++---- src/mindlessgen/qm/xtb.py | 4 ++-- test/test_config/test_config_set_attributes.py | 17 +++++++++++++++++ 4 files changed, 24 insertions(+), 7 deletions(-) diff --git a/mindlessgen.toml b/mindlessgen.toml index 538e071..8ebf861 100644 --- a/mindlessgen.toml +++ b/mindlessgen.toml @@ -42,7 +42,7 @@ element_composition = "C:2-3, H:1-2, O:1-2, N:1-*" # > CAUTION: This option is overridden by the 'element_composition' option. # > I.e., if an element is specified in 'element_composition' with an occurrence > 0, it will be added to the molecule anyway. forbidden_elements = "57-71, 81-*" -# > Set a charge for the molecule generation. Options: "none" or "" (turn it off), "int" or +# > Set a charge for the molecule generation. Options: "none" (random charge assignment), "int" or (fixed charge assignment) molecular_charge = "none" [refine] diff --git a/src/mindlessgen/molecules/miscellaneous.py b/src/mindlessgen/molecules/miscellaneous.py index d0ad736..cdf641e 100644 --- a/src/mindlessgen/molecules/miscellaneous.py +++ b/src/mindlessgen/molecules/miscellaneous.py @@ -100,11 +100,11 @@ def calculate_ligand_electrons(natoms: np.ndarray, nel: int) -> int: Calculate the number of ligand electrons in a molecule if lanthanides or actinides are within the molecule. """ f_electrons = sum( - occurence + occurrence * ( - ati - 3 + 1 - ) # subtract 3 to get the number of electrons for an Ln3+ (Ac3+) ion. - for ati, occurence in enumerate(natoms) + ati - 2 + ) # subtract 3 to get the number of electrons for an Ln3+ (Ac3+) ion and add 1 to account for the 0 based indexing. + for ati, occurrence in enumerate(natoms) if (ati in get_lanthanides() or ati in get_actinides()) ) ligand_electrons = nel - f_electrons diff --git a/src/mindlessgen/qm/xtb.py b/src/mindlessgen/qm/xtb.py index 3569d49..3fe4435 100644 --- a/src/mindlessgen/qm/xtb.py +++ b/src/mindlessgen/qm/xtb.py @@ -104,8 +104,8 @@ def singlepoint(self, molecule: Molecule, verbosity: int = 1) -> str: Perform a single-point calculation using xtb. """ super_heavy_elements = False - ati_original = molecule.ati.copy() if np.any(molecule.ati > 85): + ati_original = molecule.ati.copy() super_heavy_elements = True molecule.ati[molecule.ati > 85] -= 32 if np.any(np.isin(molecule.ati, get_lanthanides())) or np.any( @@ -153,7 +153,7 @@ def singlepoint(self, molecule: Molecule, verbosity: int = 1) -> str: molecule.uhf = uhf_original if super_heavy_elements: # Reset the atomic numbers to the original values before returning the optimized molecule. - molecule.ati = ati_original + molecule.ati = ati_original # pylint: disable=E0606 return xtb_log_out def check_gap( diff --git a/test/test_config/test_config_set_attributes.py b/test/test_config/test_config_set_attributes.py index 0af85f3..2092435 100644 --- a/test/test_config/test_config_set_attributes.py +++ b/test/test_config/test_config_set_attributes.py @@ -141,6 +141,23 @@ def test_generate_config_molecular_charge(property_name, valid_value, expected_v assert getattr(config, property_name) == expected_value +@pytest.mark.parametrize( + "property_name, invalid_value, expected_exception", + [ + ("molecular_charge", "two", ValueError), + ("molecular_charge", "1.0", ValueError), + ], +) +def test_generate_config_molecular_charge_invalid( + property_name, invalid_value, expected_exception +): + config = GenerateConfig() + + # Test invalid value + with pytest.raises(expected_exception): + setattr(config, property_name, invalid_value) + + # Tests for RefineConfig @pytest.mark.parametrize( "property_name, valid_value, invalid_value, expected_exception",