From 649b7ca5323bbd2f2fddc22a726e8f90d105add5 Mon Sep 17 00:00:00 2001 From: GES-lbabetto Date: Thu, 31 Mar 2022 10:29:58 +0200 Subject: [PATCH 1/7] added --fstrict flag to (de)protonation routines by default --- compechem/modules/crest.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/compechem/modules/crest.py b/compechem/modules/crest.py index 3b8d9dc..b256932 100644 --- a/compechem/modules/crest.py +++ b/compechem/modules/crest.py @@ -151,7 +151,7 @@ def deprotonate(mol, nproc=1, remove_tdir=True, optionals=""): mol.write_xyz("geom.xyz") os.system( - f"crest geom.xyz --alpb water --chrg {mol.charge} --uhf {mol.spin-1} --deprotonate {optionals} -T {nproc} > output.out 2>> output.err" + f"crest geom.xyz --alpb water --chrg {mol.charge} --uhf {mol.spin-1} --deprotonate --fstrict {optionals} -T {nproc} > output.out 2>> output.err" ) if os.path.exists("deprotonated.xyz"): @@ -213,7 +213,7 @@ def protonate(mol, nproc=1, remove_tdir=True, optionals=""): mol.write_xyz("geom.xyz") os.system( - f"crest geom.xyz --alpb water --chrg {mol.charge} --uhf {mol.spin-1} --protonate {optionals} -T {nproc} > output.out 2>> output.err" + f"crest geom.xyz --alpb water --chrg {mol.charge} --uhf {mol.spin-1} --protonate --fstrict {optionals} -T {nproc} > output.out 2>> output.err" ) if os.path.exists("protonated.xyz"): From 065fd4d353e9a1a629d76ec50b50e5cd1ade0cea Mon Sep 17 00:00:00 2001 From: GES-lbabetto Date: Mon, 4 Apr 2022 12:27:00 +0200 Subject: [PATCH 2/7] r2SCAN routine caused optimizations to fail with default settings --- compechem/modules/orca.py | 2 +- compechem/molecule.py | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/compechem/modules/orca.py b/compechem/modules/orca.py index 91920ea..ddc3ad1 100644 --- a/compechem/modules/orca.py +++ b/compechem/modules/orca.py @@ -410,7 +410,7 @@ def __init__(self, nproc=1, maxcore=350): super().__init__( method="r2SCAN-3c", basis_set="", - aux_basis="", + aux_basis="\n", nproc=nproc, maxcore=maxcore, solvation=True, diff --git a/compechem/molecule.py b/compechem/molecule.py index 3c4ccae..a069553 100644 --- a/compechem/molecule.py +++ b/compechem/molecule.py @@ -38,6 +38,7 @@ def __init__(self, xyz_file, charge=0, spin=1) -> None: self.spin: int = spin self.atomcount: int = None self.geometry: list = [] + self.flags: list = [] self.energies: dict = {} From 1bc474431e895c516f1fe542c5e046b24b77da3c Mon Sep 17 00:00:00 2001 From: GES-lbabetto Date: Mon, 4 Apr 2022 12:37:34 +0200 Subject: [PATCH 3/7] Found the actual bug for Orca geometry optimizations... --- compechem/modules/orca.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/compechem/modules/orca.py b/compechem/modules/orca.py index ddc3ad1..618c901 100644 --- a/compechem/modules/orca.py +++ b/compechem/modules/orca.py @@ -173,7 +173,7 @@ def opt(self, mol, charge=None, spin=None, inplace=False, remove_tdir=True): f"%pal nproc {self.nproc} end\n" f"%maxcore {self.maxcore}\n" f"! {self.method} {self.basis_set} {self.optionals}\n" - f"! RIJCOSX {self.aux_basis}" + f"! RIJCOSX {self.aux_basis}\n" ) if self.solvation is True: inp.write( @@ -410,7 +410,7 @@ def __init__(self, nproc=1, maxcore=350): super().__init__( method="r2SCAN-3c", basis_set="", - aux_basis="\n", + aux_basis="", nproc=nproc, maxcore=maxcore, solvation=True, From 8a0e1c80e5704f4856022251f41dd0ed274a9755 Mon Sep 17 00:00:00 2001 From: GES-lbabetto Date: Mon, 4 Apr 2022 14:55:28 +0200 Subject: [PATCH 4/7] fixed wrong return in Orca opt --- compechem/modules/orca.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/compechem/modules/orca.py b/compechem/modules/orca.py index 618c901..863cf2e 100644 --- a/compechem/modules/orca.py +++ b/compechem/modules/orca.py @@ -217,7 +217,8 @@ def opt(self, mol, charge=None, spin=None, inplace=False, remove_tdir=True): tools.process_output(mol, self.method, charge, spin, "opt", tdir, remove_tdir, parent_dir) - return newmol + if inplace is False: + return newmol def freq(self, mol, charge=None, spin=None, inplace=False, remove_tdir=True): """Frequency analysis (analytical frequencies). From da9fae89a93975f781abaa24390efa908e9baafd Mon Sep 17 00:00:00 2001 From: GES-lbabetto Date: Tue, 5 Apr 2022 10:59:49 +0200 Subject: [PATCH 5/7] Added flag property for Molecule object + tools for adding flags in common calculations --- compechem/ensemble.py | 6 ++- compechem/modules/algorithms.py | 9 +++- compechem/modules/crest.py | 77 ++++++++++++++++++++------------- compechem/modules/orca.py | 36 +++++++++++++-- compechem/modules/tools.py | 50 +++++++++++++++++---- compechem/modules/xtb.py | 33 ++++++++++++-- compechem/molecule.py | 8 ++-- 7 files changed, 166 insertions(+), 53 deletions(-) diff --git a/compechem/ensemble.py b/compechem/ensemble.py index 5051983..97eba30 100644 --- a/compechem/ensemble.py +++ b/compechem/ensemble.py @@ -17,7 +17,7 @@ class Ensemble: calculated at various levels of theory """ - def __init__(self, molecules_list) -> None: + def __init__(self, molecules_list: list) -> None: """ Parameters ---------- @@ -53,7 +53,9 @@ def __init__( self.electronic = electronic self.vibronic = vibronic - def boltzmann_average(self, method_el, method_vib=None, temperature=297.15): + def boltzmann_average( + self, method_el: str, method_vib: str = None, temperature: float = 297.15 + ): """Calculates the average free Gibbs energy of the ensemble (in Hartree), weighted for each molecule by its Boltzmann factor. diff --git a/compechem/modules/algorithms.py b/compechem/modules/algorithms.py index 9a36fea..f4a056a 100644 --- a/compechem/modules/algorithms.py +++ b/compechem/modules/algorithms.py @@ -1,7 +1,10 @@ from compechem.ensemble import Ensemble +from compechem.molecule import Molecule -def calculate_pka(protonated, deprotonated, method_el, method_vib=None): +def calculate_pka( + protonated: Molecule, deprotonated: Molecule, method_el: str, method_vib: str = None +): """Calculates the pKa of a molecule, given the protonated and deprotonated forms. Parameters @@ -55,7 +58,9 @@ def calculate_pka(protonated, deprotonated, method_el, method_vib=None): return pka -def calculate_potential(oxidised, reduced, method_el, method_vib=None, pH=7): +def calculate_potential( + oxidised: Molecule, reduced: Molecule, method_el: str, method_vib: str = None, pH: float = 7.0 +): """Calculates the reduction potential of a molecule, given the oxidised and reduced forms. Parameters diff --git a/compechem/modules/crest.py b/compechem/modules/crest.py index b256932..ad36611 100644 --- a/compechem/modules/crest.py +++ b/compechem/modules/crest.py @@ -1,10 +1,11 @@ import os, copy from tempfile import mkdtemp +from turtle import st from compechem.molecule import Molecule from compechem.modules import tools -def tautomer_search(mol, nproc=1, remove_tdir=True, optionals=""): +def tautomer_search(mol: Molecule, nproc: int = 1, remove_tdir: bool = True, optionals: str = ""): """Tautomer search using CREST. Parameters @@ -48,6 +49,10 @@ def tautomer_search(mol, nproc=1, remove_tdir=True, optionals=""): print( f"WARNING: cyclization change spotted for {tautomer.name}. Removing from list." ) + tools.add_flag( + mol, + f"Cyclization change occurred for {tautomer.name} during conformer search. Conformer was removed.", + ) else: tautomers.append(tautomer) @@ -61,10 +66,11 @@ def tautomer_search(mol, nproc=1, remove_tdir=True, optionals=""): tools.process_output( mol, "CREST", mol.charge, mol.spin, "tautomers", tdir, remove_tdir, parent_dir ) + tools.add_flag(mol, "No possible tautomers. Tautomer search was ignored.") return [mol] -def conformer_search(mol, nproc=1, remove_tdir=True, optionals=""): +def conformer_search(mol: Molecule, nproc: int = 1, remove_tdir: bool = True, optionals: str = ""): """Conformer search using CREST. Parameters @@ -108,21 +114,26 @@ def conformer_search(mol, nproc=1, remove_tdir=True, optionals=""): print( f"WARNING: cyclization change spotted for {conformer.name}. Removing from list." ) + tools.add_flag( + mol, + f"Cyclization change occurred for {conformer.name} during conformer search. Conformer was removed.", + ) else: conformers.append(conformer) tools.process_output( mol, "CREST", mol.charge, mol.spin, "conformers", tdir, remove_tdir, parent_dir ) - return conformers + else: print("ERROR: conformer search failed.") + tools.add_flag(mol, "Error during conformer search.") os.chdir(parent_dir) return None -def deprotonate(mol, nproc=1, remove_tdir=True, optionals=""): +def deprotonate(mol: Molecule, nproc: int = 1, remove_tdir: bool = True, optionals: str = ""): """Deprotomer search using CREST. Parameters @@ -168,6 +179,10 @@ def deprotonate(mol, nproc=1, remove_tdir=True, optionals=""): print( f"WARNING: cyclization change spotted for {deprotomer.name}. Removing from list." ) + tools.add_flag( + mol, + f"Cyclization change occurred for {deprotomer.name} during deprotomer search. Deprotomer was removed.", + ) else: deprotomers.append(deprotomer) @@ -181,10 +196,8 @@ def deprotonate(mol, nproc=1, remove_tdir=True, optionals=""): os.chdir(parent_dir) return None - return deprotomers - -def protonate(mol, nproc=1, remove_tdir=True, optionals=""): +def protonate(mol: Molecule, nproc: int = 1, remove_tdir: bool = True, optionals: str = ""): """Protomer search using CREST. Parameters @@ -230,8 +243,12 @@ def protonate(mol, nproc=1, remove_tdir=True, optionals=""): print( f"WARNING: cyclization change spotted for {protomer.name}. Removing from list." ) + tools.add_flag( + mol, + f"Cyclization change occurred for {protomer.name} during deprotomer search. Protomer was removed.", + ) else: - protomers.append(deprotomer) + protomers.append(protomer) tools.process_output( mol, "CREST", mol.charge, mol.spin, "protomers", tdir, remove_tdir, parent_dir @@ -243,19 +260,17 @@ def protonate(mol, nproc=1, remove_tdir=True, optionals=""): os.chdir(parent_dir) return None - return protomers - def qcg_grow( - solute, - solvent, - charge=None, - spin=None, - method="gfn2", - nsolv=0, - nproc=1, - optionals="", - remove_tdir=True, + solute: Molecule, + solvent: Molecule, + charge: int = None, + spin: int = None, + method: str = "gfn2", + nsolv: int = 0, + nproc: int = 1, + optionals: str = "", + remove_tdir: bool = True, ): """Quantum Cluster Growth using CREST. @@ -316,6 +331,7 @@ def qcg_grow( cluster.update_geometry("grow/cluster.xyz") except: print("ERROR: cluster growth failed.") + tools.add_flag(solute, "Cluster growth failed.") os.chdir(parent_dir) return None @@ -325,17 +341,17 @@ def qcg_grow( def qcg_ensemble( - solute, - solvent, - charge=None, - spin=None, - method="gfn2", - enslvl="gfn2", - ensemble_choice="full_ensemble", - nsolv=0, - nproc=1, - optionals="", - remove_tdir=True, + solute: Molecule, + solvent: Molecule, + charge: int = None, + spin: int = None, + method: st = "gfn2", + enslvl: str = "gfn2", + ensemble_choice: str = "full_ensemble", + nsolv: int = 0, + nproc: int = 1, + optionals: str = "", + remove_tdir: bool = True, ): """Quantum Cluster Growth + ensemble generation using CREST. @@ -406,6 +422,7 @@ def qcg_ensemble( except: print("ERROR: cluster growth failed.") + tools.add_flag(solute, "Cluster growth failed.") os.chdir(parent_dir) return None diff --git a/compechem/modules/orca.py b/compechem/modules/orca.py index 863cf2e..01f8658 100644 --- a/compechem/modules/orca.py +++ b/compechem/modules/orca.py @@ -49,7 +49,14 @@ def __init__( self.solvent = solvent self.optionals = optionals - def spe(self, mol, charge=None, spin=None, inplace=False, remove_tdir=True): + def spe( + self, + mol: Molecule, + charge: int = None, + spin: int = None, + inplace: bool = False, + remove_tdir: bool = True, + ): """Single point energy calculation. Parameters @@ -130,7 +137,14 @@ def spe(self, mol, charge=None, spin=None, inplace=False, remove_tdir=True): if inplace is False: return newmol - def opt(self, mol, charge=None, spin=None, inplace=False, remove_tdir=True): + def opt( + self, + mol: Molecule, + charge: int = None, + spin: int = None, + inplace: bool = False, + remove_tdir: bool = True, + ): """Geometry optimization + frequency analysis. Parameters @@ -220,7 +234,14 @@ def opt(self, mol, charge=None, spin=None, inplace=False, remove_tdir=True): if inplace is False: return newmol - def freq(self, mol, charge=None, spin=None, inplace=False, remove_tdir=True): + def freq( + self, + mol: Molecule, + charge: int = None, + spin: int = None, + inplace: bool = False, + remove_tdir: bool = True, + ): """Frequency analysis (analytical frequencies). Note: if the SMD solvation model is detected, defaults to numerical frequencies @@ -309,7 +330,14 @@ def freq(self, mol, charge=None, spin=None, inplace=False, remove_tdir=True): if inplace is False: return newmol - def nfreq(self, mol, charge=None, spin=None, inplace=False, remove_tdir=True): + def nfreq( + self, + mol: Molecule, + charge: int = None, + spin: int = None, + inplace: bool = False, + remove_tdir: bool = True, + ): """Frequency analysis (numerical frequencies). Parameters diff --git a/compechem/modules/tools.py b/compechem/modules/tools.py index 80b7b43..c2aee63 100644 --- a/compechem/modules/tools.py +++ b/compechem/modules/tools.py @@ -1,13 +1,14 @@ import os import shutil import pickle +from sys import flags from rdkit import Chem from compechem.molecule import Molecule from compechem.modules.orca import r2SCAN from compechem.modules.xtb import XtbInput -def generate_inchikey(molfile): +def generate_inchikey(molfile: str): """Generates InchiKey from a Mol file. Parameters @@ -25,7 +26,7 @@ def generate_inchikey(molfile): return inchikey -def generate_inchi(molfile): +def generate_inchi(molfile: str): """Generates Inchi from a Mol file. Parameters @@ -43,7 +44,21 @@ def generate_inchi(molfile): return inchi -def info(mol, print_geometry=True): +def add_flag(mol: Molecule, flag: str): + """Adds a warning flag to a Molecule object. + + Parameters + ---------- + mol : Molecule object + Molecule object to which the flag will be added + flag : str + String representing the warning which needs to be added + """ + mol.flags.append(flag) + return + + +def info(mol: Molecule, print_geometry: bool = True): """Prints information about the molecule Parameters @@ -61,6 +76,9 @@ def info(mol, print_geometry=True): print(f"\nNumber of atoms: {mol.atomcount}") print(f"Charge: {mol.charge}") print(f"Spin: {mol.spin}") + print("\n--- Warnings ---\n") + for warning in mol.flags: + print(warning) print("\n --- Energies (Eh) --- \n") for method in mol.energies: print(f"* Method: {method}") @@ -72,7 +90,7 @@ def info(mol, print_geometry=True): print(line, end="") -def dump(obj, filename=None): +def dump(obj, filename: bool = None): """Generates a pickle file containing an object Parameters @@ -93,7 +111,7 @@ def dump(obj, filename=None): pickle.dump(obj, open(filename, "wb")) -def save_ext(ext, output_dir): +def save_ext(ext: list, output_dir: str): """Saves all files matching a certain set of extensions Parameters @@ -115,7 +133,16 @@ def save_ext(ext, output_dir): shutil.copy(file, output_dir + "/" + file) -def process_output(mol, method, charge, spin, calc, tdir, remove_tdir, parent_dir): +def process_output( + mol: Molecule, + method: str, + charge: int, + spin: int, + calc: str, + tdir: str, + remove_tdir: bool, + parent_dir: str, +): """Processes the output of a calculation, copying the output files to a safe directory in the parent directory tree, and cleans the temporary directory if requested. @@ -156,7 +183,7 @@ def process_output(mol, method, charge, spin, calc, tdir, remove_tdir, parent_di os.chdir(parent_dir) -def cyclization_check(start_file, end_file): +def cyclization_check(start_file: str, end_file: str): """Checks if a cyclization has occurred (e.g., during a geometry optimization), or if a ring opening has occurred. @@ -221,7 +248,7 @@ def dissociation_check(): return False -def split_multixyz(mol, file, suffix, charge=None, spin=None): +def split_multixyz(mol: Molecule, file: str, suffix: str, charge: int = None, spin: int = None): """Splits a .xyz file containing multiple structures into individual structures. Parameters @@ -268,7 +295,12 @@ def split_multixyz(mol, file, suffix, charge=None, spin=None): def reorder_energies( - molecule_list, nproc=1, maxcore=350, method_opt=None, method_el=None, method_vib=None + molecule_list: list, + nproc: int = 1, + maxcore: int = 350, + method_opt: str = None, + method_el: str = None, + method_vib: str = None, ): """Reorders a molecule list (generated by a CREST routine such as deprotonation) at a different level of theory. diff --git a/compechem/modules/xtb.py b/compechem/modules/xtb.py index 55169f0..8327840 100644 --- a/compechem/modules/xtb.py +++ b/compechem/modules/xtb.py @@ -38,7 +38,14 @@ def __init__( self.solvent = solvent self.optionals = optionals - def spe(self, mol, charge=None, spin=None, inplace=False, remove_tdir=True): + def spe( + self, + mol: Molecule, + charge: int = None, + spin: int = None, + inplace: bool = False, + remove_tdir: bool = True, + ): """Single point energy calculation. Parameters @@ -116,7 +123,14 @@ def spe(self, mol, charge=None, spin=None, inplace=False, remove_tdir=True): if inplace is False: return newmol - def opt(self, mol, charge=None, spin=None, inplace=False, remove_tdir=True): + def opt( + self, + mol: Molecule, + charge: int = None, + spin: int = None, + inplace: bool = False, + remove_tdir: bool = True, + ): """Geometry optimization + frequency analysis. Parameters @@ -169,10 +183,16 @@ def opt(self, mol, charge=None, spin=None, inplace=False, remove_tdir=True): if tools.dissociation_check() is True: print(f"ERROR: dissociation spotted for {mol.name}.") + tools.add_flag( + mol, f"Dissociation occurred during geometry optimization with {self.method}." + ) os.chdir(parent_dir) return None elif tools.cyclization_check(f"{mol.name}.xyz", "xtbopt.xyz") is True: print(f"ERROR: cyclization change spotted for {mol.name}.") + tools.add_flag( + mol, f"Cyclization change occurred during geometry optimization with {self.method}." + ) os.chdir(parent_dir) return None else: @@ -209,7 +229,14 @@ def opt(self, mol, charge=None, spin=None, inplace=False, remove_tdir=True): if inplace is False: return newmol - def freq(self, mol, charge=None, spin=None, inplace=False, remove_tdir=True): + def freq( + self, + mol: Molecule, + charge: int = None, + spin: int = None, + inplace: bool = False, + remove_tdir: bool = True, + ): """Frequency analysis. Parameters diff --git a/compechem/molecule.py b/compechem/molecule.py index a069553..f14d080 100644 --- a/compechem/molecule.py +++ b/compechem/molecule.py @@ -19,9 +19,11 @@ class Molecule: energies : dict dictionary containing the electronic/vibronic energies of the molecule, calculated at various levels of theory + flags : list + list containing all "warning" flags which might be encountered during calculations. """ - def __init__(self, xyz_file, charge=0, spin=1) -> None: + def __init__(self, xyz_file: str, charge: int = 0, spin: int = 1) -> None: """ Parameters ---------- @@ -74,7 +76,7 @@ def __init__( def __str__(self): return f"method: {self.method}, el={self.electronic}, vib={self.vibronic}" - def write_xyz(self, xyz_file): + def write_xyz(self, xyz_file: str): """Writes the current geometry to a .xyz file. Parameters @@ -88,7 +90,7 @@ def write_xyz(self, xyz_file): for line in self.geometry: file.write(line) - def update_geometry(self, xyz_file): + def update_geometry(self, xyz_file: str): """Updates the current geometry from an external .xyz file Parameters From fd45035b972f29ccab6802ad9416b9fac56c9645 Mon Sep 17 00:00:00 2001 From: GES-lbabetto Date: Wed, 6 Apr 2022 09:33:36 +0200 Subject: [PATCH 6/7] CREST conformer search now returns original molecule in case of failed search (e.g. if a cyclization occurs) --- compechem/modules/crest.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/compechem/modules/crest.py b/compechem/modules/crest.py index ad36611..f2ce9db 100644 --- a/compechem/modules/crest.py +++ b/compechem/modules/crest.py @@ -130,7 +130,7 @@ def conformer_search(mol: Molecule, nproc: int = 1, remove_tdir: bool = True, op print("ERROR: conformer search failed.") tools.add_flag(mol, "Error during conformer search.") os.chdir(parent_dir) - return None + return [mol] def deprotonate(mol: Molecule, nproc: int = 1, remove_tdir: bool = True, optionals: str = ""): From 76d9397fd9e97cd72141138c346cbc312db83075 Mon Sep 17 00:00:00 2001 From: GES-lbabetto Date: Wed, 6 Apr 2022 11:56:39 +0200 Subject: [PATCH 7/7] version bump --- meta.yaml | 2 +- setup.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/meta.yaml b/meta.yaml index 0a74292..26fa9a7 100644 --- a/meta.yaml +++ b/meta.yaml @@ -1,5 +1,5 @@ {% set name = "GES-comp-echem" %} -{% set version = "0.1.22a" %} +{% set version = "0.1.23a" %} diff --git a/setup.py b/setup.py index bf51ca7..239912d 100644 --- a/setup.py +++ b/setup.py @@ -3,7 +3,7 @@ setuptools.setup( name="GES-comp-echem", - version="0.1.22a", + version="0.1.23a", description="", long_description="", packages=["compechem"],