Add solvation keyword to all calculation functions, move run_crest to…

… run.py and mirror run_xtb, add conformers function

conformers.py

@@ -40,25 +40,25 @@
 from shutil import rmtree, copytree
 
 from config import config, calc_dir, xtb_bin, crest_bin, config_file
-
+from run import run_crest
 
 # Disable if xtb and crest missing
 if xtb_bin is None or crest_bin is None:
     quit()
 
 
-def run_crest(command, geom_file, charge=0, multiplicity=1):
-    # Change working dir to that of geometry file to run crest correctly
-    os.chdir(geom_file.parent)
-    out_file = geom_file.with_name("output.out")
-
-    # Run in parallel
-    #os.environ["PATH"] += os.pathsep + path
+def conformers(geom_file, charge=0, multiplicity=1, solvation=None, ewin=6, hess=False):
+    unpaired_e = multiplicity - 1
+    command = [crest_bin, geom_file, "--xnam", xtb_bin, "--chrg", str(charge), "--uhf", str(unpaired_e), "--ewin", str(ewin)]
+    # Add solvation if requested
+    if solvation is not None:
+        command.append("--alpb")
+        command.append(solvation)
+    if hess:
+        command.extend(["--prop", "hess"])
 
     # Run crest from command line
-    calc = subprocess.run(command, capture_output=True, encoding="utf-8")
-    with open(out_file, "w", encoding="utf-8") as output:
-        output.write(calc.stdout)
+    calc, out_file = run_crest(command, geom_file)
 
     return geom_file.with_stem("crest_conformers")
 
@@ -79,7 +79,7 @@ def run_crest(command, geom_file, charge=0, multiplicity=1):
             "userOptions": {
                 "crest_bin": {
                     "type": "string",
-                    "label": "Location of the crest binary",
+                    "label": "Location of the CREST binary",
                     "default": str(crest_bin),
                     "order": 1.0
                 },
@@ -171,11 +171,7 @@ def run_crest(command, geom_file, charge=0, multiplicity=1):
     if args.display_name:
         print("Conformers…")
     if args.menu_path:
-        # Only show menu option if crest binary was found
-        if crest_bin is not None:
-            print("Extensions|Semi-empirical (xtb){770}")
-        else:
-            pass
+        print("Extensions|Semi-empirical (xtb){770}")
 
     if args.run_command:
         # Remove results of last calculation
@@ -198,26 +194,20 @@ def run_crest(command, geom_file, charge=0, multiplicity=1):
             with open(config_file, "w", encoding="utf-8") as config_path:
                 json.dump(config, config_path)
 
-        # First setup command to be passed
-        charge = avo_input["charge"]
-        # "Spin" passed by Avogadro is actually muliplicity so convert to n(unpaired e-)
-        spin = avo_input["spin"] - 1
-        command = [crest_bin, xyz_path, "--xnam", xtb_bin, "--chrg", str(charge), "--uhf", str(spin)]
-        if avo_input["hess"]:
-            command.extend(["--prop", "hess"])
         # crest takes energies in kcal so convert if provided in kJ (default)
         if config["energy_units"] == "kJ/mol":
             ewin_kcal = avo_input["ewin"] / 4.184
         else:
             ewin_kcal = avo_input["ewin"]
-        command.extend(["--ewin", str(ewin_kcal)])
-        if avo_input["solvent"] != "none":
-            command.extend(["--alpb", avo_input["solvent"]])
 
-        # Run calculation using formatted command and xyz file
-        conformers_path = run_crest(
-            command,
-            xyz_path
+        # Run calculation using xyz file
+        conformers_path = conformers(
+            xyz_path,
+            charge=avo_input["charge"],
+            multiplicity=avo_input["spin"],
+            solvation=avo_input["solvent"],
+            ewin=ewin_kcal,
+            hess=avo_input["hess"]
             )
 
         # Format everything appropriately for Avogadro

energy.py

@@ -37,18 +37,18 @@
 from pathlib import Path
 from shutil import rmtree
 
-from config import config, calc_dir, xtb_bin
+from config import config, calc_dir
 from run import run_xtb
 import convert
 
 
-def energy(geom_file, charge=0, multiplicity=1):
-    spin = multiplicity - 1
-    command = ["xtb", geom_file, "--chrg", str(charge), "--uhf", str(spin)]
-    # Add solvation if set globally
-    if config["solvent"] is not None:
+def energy(geom_file, charge=0, multiplicity=1, solvation=None):
+    unpaired_e = multiplicity - 1
+    command = ["xtb", geom_file, "--chrg", str(charge), "--uhf", str(unpaired_e)]
+    # Add solvation if requested
+    if solvation is not None:
         command.append("--alpb")
-        command.append(config["solvent"])
+        command.append(solvation)
     # Run xtb from command line
     calc, out_file, energy = run_xtb(command, geom_file)
     return energy
@@ -90,8 +90,9 @@ def energy(geom_file, charge=0, multiplicity=1):
         # Run calculation; returns energy as float in hartree
         energy_hartree = energy(
             xyz_path,
-            avo_input["charge"],
-            avo_input["spin"]
+            charge=avo_input["charge"],
+            multiplicity=avo_input["spin"],
+            solvation=config["solvent"],
             )
         # Convert energy to eV for Avogadro, other units for users
         energies = convert.convert_energy(energy_hartree, "hartree")

freq.py

@@ -37,18 +37,18 @@
 from pathlib import Path
 from shutil import rmtree
 
-from config import config, calc_dir, xtb_bin
+from config import config, calc_dir
 import convert
 from run import run_xtb
 
 
-def frequencies(geom_file, charge=0, multiplicity=1):
-    spin = multiplicity - 1
-    command = ["xtb", geom_file, "--hess", "--chrg", str(charge), "--uhf", str(spin)]
-    # Add solvation if set globally
-    if config["solvent"] is not None:
+def frequencies(geom_file, charge=0, multiplicity=1, solvation=None):
+    unpaired_e = multiplicity - 1
+    command = ["xtb", geom_file, "--hess", "--chrg", str(charge), "--uhf", str(unpaired_e)]
+    # Add solvation if requested
+    if solvation is not None:
         command.append("--alpb")
-        command.append(config["solvent"])
+        command.append(solvation)
     # Run xtb from command line
     calc, out_file, energy = run_xtb(command, geom_file)
 
@@ -91,8 +91,9 @@ def frequencies(geom_file, charge=0, multiplicity=1):
         # Run calculation; returns path to Gaussian file containing frequencies
         result_path = frequencies(
             xyz_path,
-            avo_input["charge"],
-            avo_input["spin"]
+            charge=avo_input["charge"],
+            multiplicity=avo_input["spin"],
+            solvation=config["solvent"]
             )
         # Currently Avogadro fails to convert the g98 file to cjson itself
         # So we have to convert output in g98 format to cjson ourselves

md.py

@@ -37,17 +37,17 @@
 from pathlib import Path
 from shutil import rmtree, copytree
 
-from config import config, calc_dir, xtb_bin
+from config import config, calc_dir
 from run import run_xtb
 
 
-def md(geom_file, input_file, charge=0, multiplicity=1):
+def md(geom_file, input_file, charge=0, multiplicity=1, solvation=None):
     spin = multiplicity - 1
     command = ["xtb", geom_file, "--input", input_file, "--omd", "--chrg", str(charge), "--uhf", str(spin)]
-    # Add solvation if set globally
-    if config["solvent"] is not None:
+    # Add solvation if requested
+    if solvation is not None:
         command.append("--alpb")
-        command.append(config["solvent"])
+        command.append(solvation)
     # Run xtb from command line
     calc, out_file, energy = run_xtb(command, geom_file)
     # Return path to trajectory file, along with energy
@@ -169,8 +169,9 @@ def md(geom_file, input_file, charge=0, multiplicity=1):
         trj_path = md(
             xyz_path,
             input_path,
-            avo_input["charge"],
-            avo_input["spin"]
+            charge=avo_input["charge"],
+            multiplicity=avo_input["spin"],
+            solvation=config["solvent"]
             )
 
         # Make sure that the calculation was successful before investing any more time

ohess.py

@@ -42,13 +42,13 @@
 from run import run_xtb
 
 
-def opt_freq(geom_file, charge=0, multiplicity=1):
+def opt_freq(geom_file, charge=0, multiplicity=1, solvation=None):
     spin = multiplicity - 1
     command = ["xtb", geom_file, "--ohess", "--chrg", str(charge), "--uhf", str(spin)]
-    # Add solvation if set globally
-    if config["solvent"] is not None:
+    # Add solvation if requested
+    if solvation is not None:
         command.append("--alpb")
-        command.append(config["solvent"])
+        command.append(solvation)
     # Run xtb from command line
     calc, out_file, energy = run_xtb(command, geom_file)
 
@@ -101,8 +101,9 @@ def opt_freq(geom_file, charge=0, multiplicity=1):
         # Run calculation; returns path to Gaussian file containing frequencies
         out_geom, out_freq, energy = opt_freq(
             xyz_path,
-            avo_input["charge"],
-            avo_input["spin"]
+            charge=avo_input["charge"],
+            multiplicity=avo_input["spin"],
+            solvation=config["solvent"]
             )
 
         # Convert frequencies

opt.py

@@ -40,13 +40,13 @@
 from run import run_xtb
 
 
-def optimize(geom_file, charge=0, multiplicity=1):
-    spin = multiplicity - 1
-    command = ["xtb", geom_file, "--opt", "--chrg", str(charge), "--uhf", str(spin)]
-    # Add solvation if set globally
-    if config["solvent"] is not None:
+def optimize(geom_file, charge=0, multiplicity=1, solvation=None):
+    unpaired_e = multiplicity - 1
+    command = ["xtb", geom_file, "--opt", "--chrg", str(charge), "--uhf", str(unpaired_e)]
+    # Add solvation if requested
+    if solvation is not None:
         command.append("--alpb")
-        command.append(config["solvent"])
+        command.append(solvation)
     # Run xtb from command line
     calc, out_file, energy = run_xtb(command, geom_file)
     # Return path to geometry file with same suffix, along with energy
@@ -88,8 +88,9 @@ def optimize(geom_file, charge=0, multiplicity=1):
         # Run calculation using xyz file
         result_path, energy = optimize(
             xyz_path,
-            avo_input["charge"],
-            avo_input["spin"]
+            charge=avo_input["charge"],
+            multiplicity=avo_input["spin"],
+            solvation=config["solvent"]
             )
         # Convert geometry
         cjson_path = convert.xyz_to_cjson(result_path)

orbitals.py

@@ -39,14 +39,14 @@
 from run import run_xtb
 
 
-def orbitals(geom_file, charge=0, multiplicity=1):
+def orbitals(geom_file, charge=0, multiplicity=1, solvation=None):
     spin = multiplicity - 1
     # Just do a single point calculation but pass molden option to get orbital printout
     command = ["xtb", geom_file, "--molden", "--chrg", str(charge), "--uhf", str(spin)]
-    # Add solvation if set globally
-    if config["solvent"] is not None:
+    # Add solvation if requested
+    if solvation is not None:
         command.append("--alpb")
-        command.append(config["solvent"])
+        command.append(solvation)
     # Run xtb from command line
     calc, out_file, energy = run_xtb(command, geom_file)
 
@@ -89,8 +89,9 @@ def orbitals(geom_file, charge=0, multiplicity=1):
         # Run calculation using xyz file
         result_path = orbitals(
             xyz_path,
-            avo_input["charge"],
-            avo_input["spin"]
+            charge=avo_input["charge"],
+            multiplicity=avo_input["spin"],
+            solvation=config["solvent"]
             )
 
         # Get molden file as string

run.py

@@ -39,11 +39,11 @@
 from pathlib import Path
 from shutil import rmtree, copytree
 
-from config import config, xtb_bin, calc_dir
+from config import config, xtb_bin, crest_bin, calc_dir
 import convert
 
 
-# All xtb commands rely on the functionality in this module
+# All xtb and crest commands rely on the functionality in this module
 # This thus effectively disables the menu command if executing would be impossible
 if xtb_bin is None:
     raise FileNotFoundError("xtb binary not found.")
@@ -55,12 +55,14 @@ def run_xtb(command, geom_file):
     # Change working dir to that of geometry file to run xtb correctly
     os.chdir(geom_file.parent)
     out_file = geom_file.with_name("output.out")
+
     # Replace xtb with xtb path
     if command[0] == "xtb":
         command[0] = xtb_bin
     # Replace <geometry file> string with geom_file
     if "<geometry_file>" in command:
         command[command.index("<geometry_file>")] = geom_file
+
     # Run xtb from command line
     calc = subprocess.run(command, capture_output=True, encoding="utf-8")
     with open(out_file, "w", encoding="utf-8") as output:
@@ -74,6 +76,30 @@ def run_xtb(command, geom_file):
     return calc, out_file, energy
 
 
+# Similarly, provide a generic function to run any crest calculation
+def run_crest(command, geom_file):
+    # Change working dir to that of geometry file to run crest correctly
+    os.chdir(geom_file.parent)
+    out_file = geom_file.with_name("output.out")
+
+    # Replace crest with crest path
+    if command[0] == "crest":
+        command[0] = crest_bin
+    # Replace <geometry file> string with geom_file
+    if "<geometry_file>" in command:
+        command[command.index("<geometry_file>")] = geom_file
+
+    # Run in parallel
+    #os.environ["PATH"] += os.pathsep + path
+    # Run crest from command line
+    calc = subprocess.run(command, capture_output=True, encoding="utf-8")
+    with open(out_file, "w", encoding="utf-8") as output:
+        output.write(calc.stdout)
+
+    # Return everything as a tuple including subprocess.CompletedProcess object
+    return calc, out_file
+
+
 def parse_energy(output_string):
     # Get final energy from output file
     # but don't convert here as not all calculation types give in same units