From eb7174c62bc8dd7b718722febe07f698ed898b00 Mon Sep 17 00:00:00 2001 From: Marcel Mueller Date: Tue, 27 Aug 2024 10:17:39 +0200 Subject: [PATCH] add Ln support by adaptive UHF recognition (#16) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * add Ln support by adaptive UHF recognition Signed-off-by: Marcel Müller * write also UHF file for the sake of consistency Signed-off-by: Marcel Müller --------- Signed-off-by: Marcel Müller --- CHANGELOG.md | 2 + src/mindlessgen/molecules/__init__.py | 10 +- .../molecules/generate_molecule.py | 63 ++------- src/mindlessgen/molecules/miscellaneous.py | 122 +++++++++++++++--- src/mindlessgen/molecules/molecule.py | 7 + src/mindlessgen/molecules/refinement.py | 5 +- src/mindlessgen/qm/xtb.py | 16 ++- test/test_generate/test_generate_molecule.py | 1 + test/test_main/test_main.py | 7 +- test/test_molecules/test_miscellaneous.py | 51 +++++++- test/test_molecules/test_refinement.py | 2 + test/test_qm/test_xtb.py | 1 + 12 files changed, 197 insertions(+), 90 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 4a5b883..5a4f6c7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,6 +10,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Comment line of `.xyz` file contains the total charge and number of unpaired electrons - Default ORCA calculation changed from r2SCAN-3c to PBE/def2-SVP - `verbosity = 3` always prints full QM output +- Adapted generation of number of unpaired electrons; thereby, support for Ln's +- Shifted group / element sorting definitions to miscellaneous ### Added - Optimization via DFT in the post-processing step diff --git a/src/mindlessgen/molecules/__init__.py b/src/mindlessgen/molecules/__init__.py index d8dcdb0..bed692b 100644 --- a/src/mindlessgen/molecules/__init__.py +++ b/src/mindlessgen/molecules/__init__.py @@ -7,17 +7,19 @@ generate_random_molecule, generate_coordinates, generate_atom_list, + check_distances, +) +from .refinement import iterative_optimization, detect_fragments +from .postprocess import postprocess_mol +from .miscellaneous import ( + set_random_charge, get_three_d_metals, get_four_d_metals, get_five_d_metals, get_lanthanides, get_alkali_metals, get_alkaline_earth_metals, - check_distances, ) -from .refinement import iterative_optimization, detect_fragments -from .postprocess import postprocess_mol -from .miscellaneous import set_random_charge __all__ = [ "Molecule", diff --git a/src/mindlessgen/molecules/generate_molecule.py b/src/mindlessgen/molecules/generate_molecule.py index b0777ac..9c2fa2d 100644 --- a/src/mindlessgen/molecules/generate_molecule.py +++ b/src/mindlessgen/molecules/generate_molecule.py @@ -6,7 +6,15 @@ import numpy as np from ..prog import GenerateConfig from .molecule import Molecule -from .miscellaneous import set_random_charge +from .miscellaneous import ( + set_random_charge, + get_alkali_metals, + get_alkaline_earth_metals, + get_three_d_metals, + get_four_d_metals, + get_five_d_metals, + get_lanthanides, +) def generate_random_molecule( @@ -29,8 +37,7 @@ def generate_random_molecule( inc_scaling_factor=config_generate.increase_scaling_factor, verbosity=verbosity, ) - mol.charge = set_random_charge(mol.ati, verbosity) - mol.uhf = 0 + mol.charge, mol.uhf = set_random_charge(mol.ati, verbosity) mol.set_name_from_formula() if verbosity > 1: @@ -359,53 +366,3 @@ def check_distances(xyz: np.ndarray, threshold: float) -> bool: if r < threshold: return False return True - - -def get_alkali_metals() -> list[int]: - """ - Get the atomic numbers of alkali metals. - """ - alkali = [2, 10, 18, 36, 54] - return alkali - - -def get_alkaline_earth_metals() -> list[int]: - """ - Get the atomic numbers of alkaline earth metals. - """ - alkaline = [3, 11, 19, 37, 55] - return alkaline - - -def get_three_d_metals() -> list[int]: - """ - Get the atomic numbers of three d metals. - """ - threedmetals = list(range(20, 30)) - - return threedmetals - - -def get_four_d_metals() -> list[int]: - """ - Get the atomic numbers of four d metals. - """ - - fourdmetals = list(range(38, 48)) - return fourdmetals - - -def get_five_d_metals() -> list[int]: - """ - Get the atomic numbers of five d metals. - """ - fivedmetals = list(range(71, 80)) - return fivedmetals - - -def get_lanthanides() -> list[int]: - """ - Get the atomic numbers of lanthanides. - """ - lanthanides = list(range(56, 71)) - return lanthanides diff --git a/src/mindlessgen/molecules/miscellaneous.py b/src/mindlessgen/molecules/miscellaneous.py index 93aa177..a735eb8 100644 --- a/src/mindlessgen/molecules/miscellaneous.py +++ b/src/mindlessgen/molecules/miscellaneous.py @@ -5,7 +5,7 @@ import numpy as np -def set_random_charge(ati: np.ndarray, verbosity: int = 1) -> int: +def set_random_charge(ati: np.ndarray, verbosity: int = 1) -> tuple[int, int]: """ Set the charge of a molecule so that unpaired electrons are avoided. """ @@ -15,22 +15,112 @@ def set_random_charge(ati: np.ndarray, verbosity: int = 1) -> int: nel += atom + 1 if verbosity > 1: print(f"Number of protons in molecule: {nel}") - iseven = False - if nel % 2 == 0: - iseven = True - # if the number of electrons is even, the charge is -2, 0, or 2 - # if the number of electrons is odd, the charge is -1, 1 - randint = np.random.rand() - if iseven: - if randint < 1 / 3: - charge = -2 - elif randint < 2 / 3: + + if np.any(np.isin(ati, get_lanthanides())): + ### Special mode for lanthanides + # -> always high spin + # -> Divide the molecule into Ln3+ ions and negative "ligands" + # -> The ligands are the remaining protons are assumed to be low spin + uhf = 0 + charge = 0 + ln_protons = 0 + for atom in ati: + if atom in get_lanthanides(): + if atom < 64: + uhf += atom - 56 + else: + uhf += 70 - atom + ln_protons += ( + atom - 3 + 1 + ) # subtract 3 to get the number of protons in the Ln3+ ion + ligand_protons = nel - ln_protons + if verbosity > 2: + print(f"Number of protons from Ln^3+ ions: {ln_protons}") + print( + f"Number of protons from ligands (assuming negative charge): {ligand_protons}" + ) + if ligand_protons % 2 == 0: charge = 0 else: - charge = 2 + charge = 1 + return charge, uhf else: - if randint < 0.5: - charge = -1 + ### Default mode + iseven = False + if nel % 2 == 0: + iseven = True + # if the number of electrons is even, the charge is -2, 0, or 2 + # if the number of electrons is odd, the charge is -1, 1 + randint = np.random.rand() + if iseven: + if randint < 1 / 3: + charge = -2 + elif randint < 2 / 3: + charge = 0 + else: + charge = 2 else: - charge = 1 - return charge + if randint < 0.5: + charge = -1 + else: + charge = 1 + uhf = 0 + return charge, uhf + + +def get_alkali_metals() -> list[int]: + """ + Get the atomic numbers of alkali metals. + """ + alkali = [2, 10, 18, 36, 54] + return alkali + + +def get_alkaline_earth_metals() -> list[int]: + """ + Get the atomic numbers of alkaline earth metals. + """ + alkaline = [3, 11, 19, 37, 55] + return alkaline + + +def get_three_d_metals() -> list[int]: + """ + Get the atomic numbers of three d metals. + """ + threedmetals = list(range(20, 30)) + + return threedmetals + + +def get_four_d_metals() -> list[int]: + """ + Get the atomic numbers of four d metals. + """ + + fourdmetals = list(range(38, 48)) + return fourdmetals + + +def get_five_d_metals() -> list[int]: + """ + Get the atomic numbers of five d metals. + """ + fivedmetals = list(range(71, 80)) + return fivedmetals + + +def get_lanthanides() -> list[int]: + """ + Get the atomic numbers of lanthanides. + """ + lanthanides = list(range(56, 71)) + return lanthanides + + +def get_actinides() -> list[int]: + """ + Get the atomic numbers of actinides. + """ + actinides = list(range(88, 103)) + return actinides diff --git a/src/mindlessgen/molecules/molecule.py b/src/mindlessgen/molecules/molecule.py index c9f5cf9..2bd16f7 100644 --- a/src/mindlessgen/molecules/molecule.py +++ b/src/mindlessgen/molecules/molecule.py @@ -340,6 +340,9 @@ def uhf(self, value: int | float): except ValueError as e: raise TypeError("Integer expected.") from e + if value < 0: + raise ValueError("Number of unpaired electrons cannot be negative.") + self._uhf = value @property @@ -483,6 +486,10 @@ def write_xyz_to_file(self, filename: str | Path | None = None): if self._charge is not None: with open(filename.with_suffix(".CHRG"), "w", encoding="utf8") as f: f.write(f"{self.charge}\n") + # if the UHF is set, write it to a '.UHF' file + if self._uhf is not None: + with open(filename.with_suffix(".UHF"), "w", encoding="utf8") as f: + f.write(f"{self.uhf}\n") def read_xyz_from_file(self, filename: str | Path): """ diff --git a/src/mindlessgen/molecules/refinement.py b/src/mindlessgen/molecules/refinement.py index 44785c5..b74b494 100644 --- a/src/mindlessgen/molecules/refinement.py +++ b/src/mindlessgen/molecules/refinement.py @@ -159,8 +159,9 @@ def detect_fragments(mol: Molecule, verbosity: int = 1) -> list[Molecule]: for atom in fragment_molecule.ati: fragment_molecule.atlist[atom] += 1 # Update the charge of the fragment molecule - fragment_molecule.charge = set_random_charge(fragment_molecule.ati, verbosity) - fragment_molecule.uhf = 0 + fragment_molecule.charge, fragment_molecule.uhf = set_random_charge( + fragment_molecule.ati, verbosity + ) fragment_molecule.set_name_from_formula() if verbosity > 1: print(f"Fragment molecule: {fragment_molecule}") diff --git a/src/mindlessgen/qm/xtb.py b/src/mindlessgen/qm/xtb.py index f1e6c5a..f42cc63 100644 --- a/src/mindlessgen/qm/xtb.py +++ b/src/mindlessgen/qm/xtb.py @@ -47,11 +47,14 @@ def optimize( "--opt", "--gfn", "2", - "--chrg", - str(molecule.charge), ] + if molecule.charge != 0: + arguments += ["--chrg", str(molecule.charge)] + if molecule.uhf != 0: + arguments += ["--uhf", str(molecule.uhf)] if max_cycles is not None: arguments += ["--cycles", str(max_cycles)] + if verbosity > 2: print(f"Running command: {' '.join(arguments)}") @@ -87,10 +90,13 @@ def singlepoint(self, molecule: Molecule, verbosity: int = 1) -> str: "molecule.xyz", "--gfn", "2", - "--chrg", - str(molecule.charge), ] - if verbosity > 1: + if molecule.charge != 0: + arguments += ["--chrg", str(molecule.charge)] + if molecule.uhf != 0: + arguments += ["--uhf", str(molecule.uhf)] + + if verbosity > 2: print(f"Running command: {' '.join(arguments)}") xtb_log_out, xtb_log_err, return_code = self._run( diff --git a/test/test_generate/test_generate_molecule.py b/test/test_generate/test_generate_molecule.py index 2116239..929041d 100644 --- a/test/test_generate/test_generate_molecule.py +++ b/test/test_generate/test_generate_molecule.py @@ -181,6 +181,7 @@ def test_generate_molecule() -> None: """ # create a ConfigManager object with verbosity set to 0 config = ConfigManager() + config.generate.forbidden_elements = "57-71" config.general.verbosity = 0 mol = generate_random_molecule(config.generate, config.general.verbosity) diff --git a/test/test_main/test_main.py b/test/test_main/test_main.py index d6faf15..98d6d31 100644 --- a/test/test_main/test_main.py +++ b/test/test_main/test_main.py @@ -8,10 +8,11 @@ @pytest.mark.optional def test_generator(): config = ConfigManager() - config.general.engine = "xtb" + config.refine.engine = "xtb" config.general.max_cycles = 10000 - config.general.parallel = 8 - config.general.verbosity = 0 + config.general.parallel = 4 + config.general.verbosity = -1 + config.general.postprocess = False molecules, exitcode = generator(config) assert exitcode == 0 diff --git a/test/test_molecules/test_miscellaneous.py b/test/test_molecules/test_miscellaneous.py index e2bdd94..68b68e5 100644 --- a/test/test_molecules/test_miscellaneous.py +++ b/test/test_molecules/test_miscellaneous.py @@ -9,15 +9,52 @@ # CAUTION: We use 0-based indexing for atoms and molecules! @pytest.mark.parametrize( - "atom_types, expected_charge", + "atom_types, expected_charges, expected_uhf", [ - (np.array([5, 7, 0]), [-1, 1]), - (np.array([5, 7, 0, 0]), [-2, 0, 2]), + # Standard mode tests (no lanthanides) + (np.array([4, 6, 0]), [-1, 1], 0), # B, N, H (odd number of electrons) + ( + np.array([4, 6, 0, 0]), + [-2, 0, 2], + 0, + ), # B, N, H, H (even number of electrons) + (np.array([5, 7, 0]), [-1, 1], 0), # C, O, H (odd number of electrons) + ( + np.array([5, 7, 0, 0]), + [-2, 0, 2], + 0, + ), # C, O, H, H (even number of electrons) + # Lanthanide mode tests + (np.array([56, 7, 0]), [0], 0), # La (57), O, H -> Lanthanide case, UHF = 0 + ( + np.array([59, 7, 7, 0]), + [0], + 3, + ), # Nd (60), O, O, H -> Lanthanide case, UHF = 3 + ( + np.array([63, 6, 6]), + [1], + 7, + ), # Gd (64), C, C -> Lanthanide case, UHF = 7, CHRG = 1 + (np.array([57, 7, 0]), [0], 1), # Ce (58), O, H -> Lanthanide case, UHF = 1 + (np.array([69, 7, 0]), [0], 1), # Ce (58), O, H -> Lanthanide case, UHF = 1 + ], + ids=[ + "B-N-H (standard, odd)", + "B-N-H-H (standard, even)", + "C-O-H (standard, odd)", + "C-O-H-H (standard, even)", + "La-O-H (lanthanide)", + "Nd-O-O-H (lanthanide)", + "Gd-C-N (lanthanide)", + "Ce-O-H (lanthanide)", + "Yb-O-H (lanthanide)", ], - ids=["odd", "even"], ) -def test_set_random_charge(atom_types, expected_charge): +def test_set_random_charge(atom_types, expected_charges, expected_uhf): """ - Test the set_random_charge function. + Test the set_random_charge function for both standard and lanthanide modes. """ - assert set_random_charge(atom_types) in expected_charge + charge, uhf = set_random_charge(atom_types) + 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 6f0250f..fba2564 100644 --- a/test/test_molecules/test_refinement.py +++ b/test/test_molecules/test_refinement.py @@ -24,6 +24,8 @@ def mol_H6O2B2Ne2I1Os1Tl1() -> Molecule: testsdir = Path(__file__).resolve().parents[1] xyz_file = testsdir / "fixtures/H6O2B2Ne2I1Os1Tl1.xyz" mol.read_xyz_from_file(xyz_file) + mol.charge = 0 + mol.uhf = 0 return mol diff --git a/test/test_qm/test_xtb.py b/test/test_qm/test_xtb.py index 3808c04..bffe308 100644 --- a/test/test_qm/test_xtb.py +++ b/test/test_qm/test_xtb.py @@ -28,6 +28,7 @@ def test_xtb_optimize_xtb(coordinates_ethanol: np.ndarray) -> None: mol.xyz = coordinates_ethanol mol.ati = np.array([7, 5, 5, 0, 0, 0, 0, 0, 0]) mol.charge = 0 + mol.uhf = 0 mol.num_atoms = 9 optimized_molecule = xtb.optimize(mol)