Merge pull request #32 from GES-lbabetto/main

Various improvements and bug fixes

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.
         

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

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
+        return [mol]
 
 
-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
@@ -151,7 +162,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"):
@@ -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
@@ -213,7 +226,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"):
@@ -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
 

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
@@ -173,7 +187,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(
@@ -217,9 +231,17 @@ 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):
+    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
@@ -308,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