Merge pull request #124 from jvalegre/jv_branch

v1.4.4
jvalegre · Feb 13, 2023 · 06776ad · 06776ad
2 parents 08812ff + f13b24a
commit 06776ad
Show file tree

Hide file tree

Showing 17 changed files with 1,560 additions and 87 deletions.
diff --git a/Example_workflows/QCORR_processing_QM_outputs/QCORR_5/success/json_files/CH4_2.json b/Example_workflows/QCORR_processing_QM_outputs/QCORR_5/success/json_files/CH4_2.json
@@ -25,7 +25,7 @@
                 "homo": -12.262538558932,
                 "gap": 17.6239977553335
             },
-            "total": -1095.6654020290605
+            "total": -1095.6754020290605
         },
         "orbitals": {
             "homos": [4],

diff --git a/Example_workflows/QCORR_processing_QM_outputs/QCORR_6/CH4_dup.log b/Example_workflows/QCORR_processing_QM_outputs/QCORR_6/CH4_dup.log
diff --git a/Example_workflows/QCORR_processing_QM_outputs/QCORR_6/success/json_files/CH4.json b/Example_workflows/QCORR_processing_QM_outputs/QCORR_6/success/json_files/CH4.json
@@ -0,0 +1,104 @@
+{
+    "chemical json": 0,
+    "name": "CH4",
+    "smiles": "C\tCH4.log\r\n",
+    "inchi": "InChI=1S/CH4/h1H4\n",
+    "inchikey": "VNWKTOKETHGBQD-UHFFFAOYSA-N\n",
+    "formula": "CH4",
+    "properties": {
+        "partial charges": {
+            "mulliken": [-0.786231, 0.196384, 0.196537, 0.196559, 0.19675],
+            "mulliken_sum": [0.0, 0.0, 0.0, 0.0, 0.0],
+            "apt": [-0.105789, 0.026726, 0.026442, 0.026078, 0.026543],
+            "apt_sum": [0.0, 0.0, 0.0, 0.0, 0.0]
+        },
+        "charge": 0,
+        "enthalpy": -40.215413,
+        "entropy": 7.870870367266035e-05,
+        "energy": {
+            "free energy": -40.23888,
+            "alpha": {
+                "homo": -12.262538558932,
+                "gap": 17.6239977553335
+            },
+            "beta": {
+                "homo": -12.262538558932,
+                "gap": 17.6239977553335
+            },
+            "total": -1095.6654020290605
+        },
+        "orbitals": {
+            "homos": [4],
+            "energies": [ [-284.6895921008575, -21.216172695784003, -12.275600023756, -12.26689238054, -12.262538558932, 5.3614591964015, 6.673864297363, 6.677129663569, 6.681483485177, 20.602011735205497, 20.6052771014115, 20.607726126066, 29.871570052488003, 30.838662677165, 30.846553978829505, 30.8634250375605, 46.7146010568865] ],
+            "molecular orbital symmetry": [
+                [ "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A" ]
+            ],
+            "basis number": 17,
+            "MO number": 17
+        },
+        "total dipole moment": 0.0008306623862918075,
+        "multiplicity": 1,
+        "number of atoms": 5,
+        "rotational": {
+            "symmetry point group": "C1",
+            "symmetry number": 1,
+            "rotational constants": [157.83492, 157.70617, 157.67138],
+            "rotational temperatures": [7.57487, 7.56869, 7.56702]
+        },
+        "pressure": 1.0,
+        "temperature": 298.15,
+        "energies": {
+            "zero-point correction": 0.045747
+        },
+        "molecular mass": 16.04246
+    },
+    "atoms": {
+        "coords": {
+            "3d": [ -0.000371, 1.5e-05, 6.6e-05, 0.457036, -0.463921, -0.87651, 0.336795, -0.516385, 0.900906, 0.29552, 1.049842, 0.050789, -1.087125, -0.069629, -0.075582 ]
+        },
+        "elements": {
+            "number": [6, 1, 1, 1, 1],
+            "atom count": 5,
+            "heavy atom count": 1
+        },
+        "core electrons": [0, 0, 0, 0, 0]
+    },
+    "metadata": {
+        "basis set": "3-21G",
+        "ground or transition state": "ground_state",
+        "dispersion model": "none",
+        "functional": "M062X",
+        "grid type": "sg1",
+        "keywords line": "opt freq 3-21g m062x",
+        "QM program": "Gaussian 09, Revision A.02",
+        "run date": "16-Jan-2022",
+        "solvation": "gas_phase"
+    },
+    "optimization": {
+        "geometric targets": [0.00045, 0.0003, 0.0018, 0.0012],
+        "geometric values": [[0.014371, 0.004905, 0.047373, 0.02206], [0.001512, 0.000647, 0.008022, 0.004003], [0.000323, 0.000117, 0.001587, 0.000724], [0.000323, 0.000117, 0.001519, 0.000674]],
+        "done": true,
+        "status": [1, 0, 4, 5],
+        "scf": {
+            "scf energies": [-1095.6422756470668, -1095.6647757617548, -1095.6654020290605, -1095.6654020290605],
+            "targets": [[1e-08, 1e-06, 1e-06], [1e-08, 1e-06, 1e-06], [1e-08, 1e-06, 1e-06], [1e-08, 1e-06, 1e-06]]
+        }
+    },
+    "vibrations": {
+        "displacement": [[[0.08, -0.07, 0.07], [-0.36, 0.32, -0.36], [-0.17, 0.04, 0.18], [-0.43, 0.11, -0.36], [0.05, 0.35, -0.32]], [[0.07, 0.1, 0.03], [-0.4, -0.4, 0.04], [-0.34, -0.49, -0.15], [-0.15, 0.13, -0.1], [0.09, -0.47, -0.09]], [[-0.07, 0.02, 0.1], [0.07, 0.17, 0.05], [0.43, -0.24, -0.26], [0.36, -0.07, -0.47], [0.0, -0.14, -0.51]], [[0.0, 0.0, 0.0], [0.35, -0.2, 0.29], [0.14, -0.39, -0.27], [-0.45, 0.13, -0.19], [-0.04, 0.47, 0.17]], [[0.0, 0.0, 0.0], [0.3, 0.41, -0.06], [-0.45, -0.19, 0.06], [0.19, -0.03, -0.46], [-0.02, -0.17, 0.47]], [[-0.01, 0.01, 0.0], [0.23, -0.23, -0.44], [0.16, -0.25, 0.45], [0.12, 0.45, 0.02], [-0.45, -0.03, -0.03]], [[0.0, -0.01, 0.09], [0.28, -0.29, -0.51], [-0.24, 0.36, -0.61], [0.03, 0.1, 0.03], [-0.02, -0.01, 0.02]], [[0.04, -0.08, -0.01], [-0.22, 0.21, 0.43], [-0.1, 0.15, -0.29], [0.2, 0.66, 0.03], [-0.34, -0.04, -0.03]], [[0.08, 0.04, 0.01], [-0.01, 0.05, 0.07], [-0.02, 0.07, -0.11], [-0.13, -0.52, -0.03], [-0.82, -0.04, -0.06]]],
+        "frequencies": [1415.554, 1418.958, 1422.4103, 1638.5106, 1641.802, 3054.0303, 3153.6387, 3162.8459, 3172.8647],
+        "force constants": [1.3958, 1.4023, 1.4088, 1.5942, 1.6007, 5.5416, 6.439, 6.4738, 6.5184],
+        "intensities": {
+            "IR": [27.4995, 27.3074, 27.4236, 0.0057, 0.0155, 0.0946, 14.4084, 14.6105, 14.4819]
+        },
+        "reduced masses": [1.1822, 1.1821, 1.1818, 1.0079, 1.0079, 1.0084, 1.0989, 1.0984, 1.099],
+        "vibration symmetry": [ "A", "A", "A", "A", "A", "A", "A", "A", "A" ]
+    },
+    "bonds": {
+        "connections": {
+            "index": [ 1, 0, 4, 0, 0, 3, 0, 2 ]
+        },
+        "order": [ 1, 1, 1, 1 ]
+    },
+    "diagram": "<?xml version=\"1.0\"?>\n<svg version=\"1.1\" id=\"topsvg\"\nxmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\"\nxmlns:cml=\"http://www.xml-cml.org/schema\" x=\"0\" y=\"0\" width=\"200px\" height=\"200px\" viewBox=\"0 0 100 100\">\n<title>CH4.log - Open Babel Depiction</title>\n<rect x=\"0\" y=\"0\" width=\"100\" height=\"100\" fill=\"white\"/>\n<g transform=\"translate(0,0)\">\n<svg width=\"100\" height=\"100\" x=\"0\" y=\"0\" viewBox=\"0 0 154.641 154.641\"\nfont-family=\"sans-serif\" stroke=\"rgb(0,0,0)\" stroke-width=\"2\"  stroke-linecap=\"round\">\n<line x1=\"103.4\" y1=\"93.5\" x2=\"80.0\" y2=\"80.0\" opacity=\"1.0\" stroke=\"rgb(0,0,0)\"  stroke-width=\"2.0\"/>\n<line x1=\"53.0\" y1=\"80.0\" x2=\"80.0\" y2=\"80.0\" opacity=\"1.0\" stroke=\"rgb(0,0,0)\"  stroke-width=\"2.0\"/>\n<line x1=\"80.0\" y1=\"80.0\" x2=\"66.5\" y2=\"103.4\" opacity=\"1.0\" stroke=\"rgb(0,0,0)\"  stroke-width=\"2.0\"/>\n<line x1=\"80.0\" y1=\"80.0\" x2=\"80.0\" y2=\"53.0\" opacity=\"1.0\" stroke=\"rgb(0,0,0)\"  stroke-width=\"2.0\"/>\n<text x=\"108.641016\" y=\"108.000000\" fill=\"rgb(191,191,191)\" stroke-width=\"0\" font-weight=\"bold\" font-size=\"16\" >H</text>\n<text x=\"74.000000\" y=\"48.000000\" fill=\"rgb(191,191,191)\" stroke-width=\"0\" font-weight=\"bold\" font-size=\"16\" >H</text>\n<text x=\"54.000000\" y=\"122.641016\" fill=\"rgb(191,191,191)\" stroke-width=\"0\" font-weight=\"bold\" font-size=\"16\" >H</text>\n<text x=\"34.000000\" y=\"88.000000\" fill=\"rgb(191,191,191)\" stroke-width=\"0\" font-weight=\"bold\" font-size=\"16\" >H</text>\n</svg>\n</g>\n<text font-size=\"18.000000\" fill =\"black\" font-family=\"sans-serif\"\nx=\"10.000000\" y=\"20.000000\" >CH4.log</text>\n</svg>\n"
+}
diff --git a/README.md b/README.md
@@ -25,8 +25,9 @@ Don't miss out the latest hands-on tutorials from our [YouTube channel](https://
 
 ## Installation
 Check our [AQME installation in 2 mins](https://youtu.be/VeaBzqIZHbo) video for a quick installation guide. In a nutshell, AQME and its dependencies are installed as follows:  
-1. Using conda-forge: `conda install -c conda-forge aqme`  
-2. Using pip: `pip install aqme`  
+1. Install using conda-forge: `conda install -c conda-forge aqme`  
+2. Install using pip: `pip install aqme`  
+3. Update to the latest version: `pip install aqme --upgrade`  
 
 Extra requirements if xTB or CREST are used (MacOS and Linux only):  
   * xTB: `conda install -y -c conda-forge xtb`  

diff --git a/README.rst b/README.rst
@@ -92,18 +92,24 @@ Check our `AQME installation in 2 mins <https://youtu.be/VeaBzqIZHbo>`_ video
 for a quick installation guide. In a nutshell, AQME and its dependencies are 
 installed as follows:
 
-1. Using conda-forge  
+1. Install using conda-forge  
 
 .. code-block:: shell 
    
    conda install -c conda-forge aqme
 
-2. Using pip:  
+2. Install using pip:  
 
 .. code-block:: shell
 
    pip install aqme
 
+3. Update to the latest version:  
+
+.. code-block:: shell
+
+   pip install aqme --upgrade
+
 Installation of the extra requirements
 ++++++++++++++++++++++++++++++++++++++
 

diff --git a/aqme/cmin.py b/aqme/cmin.py
@@ -135,7 +135,7 @@ def __init__(self, **kwargs):
             if self.args.program.lower() not in ["xtb", "ani"]:
                 cmin_program = False
         if not cmin_program:
-            self.args.log.write("\nx  Program not supported for CMIN refinement! Specify: program='xtb' (or 'ani')")
+            self.args.log.write('\nx  Program not supported for CMIN refinement! Specify: program="xtb" (or "ani")')
             self.args.log.finalize()
             sys.exit()
 

diff --git a/aqme/csearch/base.py b/aqme/csearch/base.py
@@ -227,27 +227,27 @@ def __init__(self, **kwargs):
         # load default and user-specified variables
         self.args = load_variables(kwargs, "csearch")
 
-        if self.args.program.lower() == "crest":
-            try:
-                subprocess.run(
-                    ["xtb", "-h"], stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL
-                )
-            except FileNotFoundError:
-                self.args.log.write("x  xTB is not installed (CREST cannot be used)! You can install the program with 'conda install -c conda-forge xtb'")
-                self.args.log.finalize()
-                sys.exit()
-
         csearch_program = True
         if self.args.program is None:
             csearch_program = False
         if csearch_program:
             if self.args.program.lower() not in ["rdkit", "summ", "fullmonte", "crest"]:
                 csearch_program = False
         if not csearch_program:
-            self.args.log.write("\nx  Program not supported for CSEARCH conformer generation! Specify: program='rdkit' (or summ, fullmonte, crest)")
+            self.args.log.write('\nx  Program not specified or not supported for CSEARCH! Specify: program="rdkit" (or "crest", "summ", "fullmonte")')
             self.args.log.finalize()
             sys.exit()
 
+        if self.args.program.lower() == "crest":
+            try:
+                subprocess.run(
+                    ["xtb", "-h"], stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL
+                )
+            except FileNotFoundError:
+                self.args.log.write("x  xTB is not installed (CREST cannot be used)! You can install the program with 'conda install -c conda-forge xtb'")
+                self.args.log.finalize()
+                sys.exit()
+
         if self.args.smi is None and self.args.input == "":
             self.args.log.write("\nx  Program requires either a SMILES or an input file to proceed! Please look up acceptable file formats. Specify: smi='CCC' (or input='filename.csv')")
             self.args.log.finalize()
@@ -450,12 +450,6 @@ def compute_confs(
             # check if the optimization is constrained
             complex_ts = check_constraints(self)
 
-        # this part allows to use charges and multiplicity from COM/GJF inputs
-        if self.args.charge is None:
-            self.args.charge = charge
-        if self.args.mult is None:
-            self.args.mult = mult
-
         if self.args.destination is None:
             self.csearch_folder = Path(self.args.initial_dir).joinpath(
                 f"CSEARCH"
@@ -544,6 +538,8 @@ def compute_confs(
                             constraints_angle,
                             constraints_dihedral,
                             complex_ts,
+                            charge,
+                            mult,
                             coord_map,
                             alg_map,
                             template,
@@ -562,7 +558,9 @@ def compute_confs(
                 constraints_dist,
                 constraints_angle,
                 constraints_dihedral,
-                complex_ts
+                complex_ts,
+                charge,
+                mult,
             )
 
         # Updates the dataframe with infromation about conformer generation
@@ -578,6 +576,8 @@ def conformer_generation(
         constraints_angle,
         constraints_dihedral,
         complex_ts,
+        charge,
+        mult,
         coord_Map=None,
         alg_Map=None,
         mol_template=None,
@@ -609,17 +609,13 @@ def conformer_generation(
         # Set charge and multiplicity
         metal_found = False
         # user can overwrite charge and mult with the corresponding arguments
-        if self.args.charge is not None:
-            charge = self.args.charge
-        else:
+        if charge is None:
             if not len(self.args.metal_atoms) >= 1:
                 charge = Chem.GetFormalCharge(mol)
             else:
                 charge, metal_found = rules_get_charge(mol, self.args)
 
-        if self.args.mult is not None:
-            mult = self.args.mult
-        else:
+        if mult is None:
             mult = Descriptors.NumRadicalElectrons(mol) + 1
             if metal_found:
                 # since RDKit gets the multiplicity of the metal with valence 0, the real multiplicity
@@ -757,29 +753,38 @@ def summ_search(
                     )
 
         if self.args.program.lower() in ['crest']:
+            stop_xtb_opt = False
             if not complex_ts:
                 # mol_crest is the RDKit-optimized mol object
-                rdmolfiles.MolToXYZFile(mol_crest, name + "_crest.xyz")
+                if mol_crest is not None:
+                    rdmolfiles.MolToXYZFile(mol_crest, name + "_crest.xyz")
+                else:
+                    stop_xtb_opt = True
+                    status = -1
             else:
                 # mol is the raw mol object (no optimization with RDKit to avoid problems when using
                 # noncovalent complexes and TSs)
-                rdmolfiles.MolToXYZFile(mol, name + "_crest.xyz")
-
-            status = xtb_opt_main(
-                    f"{name}_{self.args.program.lower()}",
-                    dup_data,
-                    dup_data_idx,
-                    self,
-                    charge,
-                    mult,
-                    constraints_atoms,
-                    constraints_dist,
-                    constraints_angle,
-                    constraints_dihedral,
-                    'crest',
-                    complex_ts=complex_ts,
-                    mol=mol # this is necessary for CREST calculations with constraints
-                )
+                if mol is not None:
+                    rdmolfiles.MolToXYZFile(mol, name + "_crest.xyz")
+                else:
+                    stop_xtb_opt = True
+                    status = -1
+            if not stop_xtb_opt:
+                status = xtb_opt_main(
+                        f"{name}_{self.args.program.lower()}",
+                        dup_data,
+                        dup_data_idx,
+                        self,
+                        charge,
+                        mult,
+                        constraints_atoms,
+                        constraints_dist,
+                        constraints_angle,
+                        constraints_dihedral,
+                        'crest',
+                        complex_ts=complex_ts,
+                        mol=mol # this is necessary for CREST calculations with constraints
+                    )
 
         return status
 
@@ -1193,8 +1198,13 @@ def rdkit_to_sdf(
         Conversion from RDKit to SDF
         """
 
-        Chem.SanitizeMol(mol)
-        mol = Chem.AddHs(mol)
+        try:
+            Chem.SanitizeMol(mol)
+            mol = Chem.AddHs(mol)
+        except Chem.AtomValenceException: # this happens sometimes with complex metals when substituting the metal with an I atom
+            self.args.log.write(f'\nx  The species provided could not be converted into a mol object wth RDKit. It normally happens with tricky metal complexes and might be fixed with different structures (i.e., changing a single bond + positive charge with a double bond).')
+            return -1, None, None, None
+
         mol.SetProp("_Name", name)
 
         # detects and applies auto-detection of initial number of conformers

diff --git a/aqme/csearch/utils.py b/aqme/csearch/utils.py
@@ -226,11 +226,28 @@ def generate_mol_from_csv(args, csv_smiles, index):
         constraints_angle,
         constraints_dihedral,
     )
+
+    charge_found, mult_found = False, False
+    if args.charge is None or args.mult is None:
+        for csv_column in csv_smiles.columns:
+            if str(csv_smiles.loc[index, csv_column]) != 'nan':
+                if csv_column.lower() in ['charge','chrg'] and args.charge is None:
+                    charge = int(csv_smiles.loc[index, csv_column])
+                    charge_found = True
+                elif csv_column.lower() in ['mult','multiplicity'] and args.mult is None:
+                    mult = int(csv_smiles.loc[index, csv_column])
+                    mult_found = True
+
+    if not charge_found:
+        charge = args.charge
+    if not mult_found:
+        mult = args.mult
+
     obj = (
         smiles,
         name,
-        args.charge,
-        args.mult,
+        charge,
+        mult,
         constraints_atoms,
         constraints_dist,
         constraints_angle,
@@ -484,8 +501,8 @@ def smi_to_mol(
         try:
             mol = Chem.MolFromSmiles(smi[0], params)
         except Chem.AtomValenceException:
-            log.write(f"\nx  The SMILES string provided ( {smi[0]} ) contains errors. For example, N atoms from ligands of metal complexes should be N+ since they're drawn with four bonds in ChemDraw, same for O atoms in carbonyl ligands, etc.\n")
-            sys.exit()
+            log.write(f"\nx  The SMILES string provided ( {smi[0]} ) contains errors or the molecule needs to be drawn in a different way. For example, N atoms from ligands of metal complexes should be N+ since they're drawn with four bonds in ChemDraw, same for O atoms in carbonyl ligands, etc.\n")
+            mol = None
 
     return (
         mol,