diff --git a/avo_xtb/about.py b/avo_xtb/about.py index 2bb1f00..e0ef295 100644 --- a/avo_xtb/about.py +++ b/avo_xtb/about.py @@ -4,6 +4,7 @@ import argparse import json import logging +import re import subprocess import sys @@ -13,6 +14,30 @@ logger = logging.getLogger(__name__) +# Get xtb and crest versions +# Regex to match version numbers (e.g., 6.7.1, 2.12, 6.4.0) +pattern = r"\bversion\s+(\d+\.\d+\.\d+|\d+\.\d+)\b" + +if easyxtb.XTB_BIN: + xtb_version_info = subprocess.run( + [str(easyxtb.XTB_BIN), "--version"], + encoding="utf-8", + capture_output=True, + ).stdout + XTB_VERSION = re.findall(pattern, xtb_version_info, re.IGNORECASE)[0] +else: + XTB_VERSION = None +if easyxtb.CREST_BIN: + crest_version_info = subprocess.run( + [str(easyxtb.CREST_BIN), "--version"], + encoding="utf-8", + capture_output=True, + ).stdout + CREST_VERSION = re.findall(pattern, crest_version_info, re.IGNORECASE)[0] +else: + CREST_VERSION = None + + if __name__ == "__main__": parser = argparse.ArgumentParser() parser.add_argument("--debug", action="store_true") @@ -37,30 +62,22 @@ avo_input = json.loads(sys.stdin.read()) output = avo_input.copy() - if easyxtb.XTB_BIN: - xtb_version = subprocess.run( - [str(easyxtb.XTB_BIN), "--version"], - encoding="utf-8", - capture_output=True, - ).stdout.splitlines()[-2].strip() + if XTB_VERSION: + xtb_version_msg = XTB_VERSION else: - xtb_version = "No xtb binary found" - if easyxtb.CREST_BIN: - crest_version = subprocess.run( - [str(easyxtb.CREST_BIN), "--version"], - encoding="utf-8", - capture_output=True, - ).stdout.splitlines()[-4].strip() + xtb_version_msg = "No xtb binary found" + if CREST_VERSION: + crest_version_msg = CREST_VERSION else: - crest_version = "No CREST binary found" + crest_version_msg = "No CREST binary found" # Do nothing to data other than add message with version and path info output["message"] = ( "avo_xtb plugin\n" + f"easyxtb version: {easyxtb.configuration.easyxtb_VERSION}\n" - + f"xtb version: {xtb_version}\n" + + f"xtb version: {xtb_version_msg}\n" + f"xtb path: {easyxtb.XTB_BIN}\n" - + f"CREST version: {crest_version}\n" + + f"CREST version: {crest_version_msg}\n" + f"CREST path: {easyxtb.CREST_BIN}" ) diff --git a/avo_xtb/config.py b/avo_xtb/config.py index 7380261..34b4752 100644 --- a/avo_xtb/config.py +++ b/avo_xtb/config.py @@ -33,25 +33,31 @@ "default": str(easyxtb.XTB_BIN), "order": 1.0, }, + "crest_bin": { + "type": "string", + "label": "Location of the CREST binary", + "default": str(easyxtb.CREST_BIN), + "order": 2.0, + }, "user_dir": { "type": "string", "label": "Run calculations (in subfolder) in", "default": str(easyxtb.CALC_DIR), - "order": 2.0, + "order": 3.0, }, "n_proc": { "type": "integer", "label": "Parallel processes to run", "minimum": 1, "default": 1, - "order": 3.0 + "order": 4.0 }, "energy_units": { "type": "stringList", "label": "Preferred energy units", "values": ["kJ/mol", "kcal/mol"], "default": 0, - "order": 4.0, + "order": 5.0, }, "solvent": { "type": "stringList", @@ -84,14 +90,14 @@ "water", ], "default": 0, - "order": 5.0, + "order": 6.0, }, "method": { "type": "stringList", - "label": "Method (xtb only)", + "label": "Method", "values": methods, "default": methods[-1], - "order": 6.0, + "order": 7.0, }, "opt_lvl": { "type": "stringList", @@ -107,7 +113,7 @@ "extreme", ], "default": 4, - "order": 7.0, + "order": 8.0, }, "warning": { "type": "text", @@ -146,9 +152,18 @@ easyxtb.config["calc_dir"] = str(easyxtb.CALC_DIR) # Save change to xtb_bin if there has been one - if avo_input["xtb_bin"] != str(easyxtb.XTB_BIN): + if avo_input["xtb_bin"] in ["none", ""]: + pass + elif avo_input["xtb_bin"] != str(easyxtb.XTB_BIN): easyxtb.XTB_BIN = Path(avo_input["xtb_bin"]) easyxtb.config["xtb_bin"] = str(easyxtb.XTB_BIN) + + # Save change to crest_bin if there has been one + if avo_input["crest_bin"] in ["none", ""]: + pass + elif Path(avo_input["crest_bin"]) != easyxtb.CREST_BIN: + easyxtb.CREST_BIN = Path(avo_input["crest_bin"]) + easyxtb.config["crest_bin"] = str(easyxtb.CREST_BIN) # Update number of threads easyxtb.config["n_proc"] = avo_input["n_proc"] @@ -156,14 +171,11 @@ # Update energy units easyxtb.config["energy_units"] = avo_input["energy_units"] - # Convert "none" string to Python None + # Update solvent if avo_input["solvent"] == "none": - solvent_selected = None + easyxtb.config["solvent"] = None else: - solvent_selected = avo_input["solvent"] - - # Update solvent - easyxtb.config["solvent"] = solvent_selected + easyxtb.config["solvent"] = avo_input["solvent"] # Update method easyxtb.config["method"] = methods.index(avo_input["method"]) diff --git a/avo_xtb/conformers.py b/avo_xtb/conformers.py index 2eafbd7..bc2796b 100644 --- a/avo_xtb/conformers.py +++ b/avo_xtb/conformers.py @@ -31,74 +31,7 @@ if args.print_options: options = { - "inputMoleculeFormat": "xyz", "userOptions": { - "crest_bin": { - "type": "string", - "label": "Location of the CREST binary", - "default": str(easyxtb.CREST_BIN), - "order": 1.0, - }, - "save_dir": { - "type": "string", - "label": "Save results in", - "default": str(easyxtb.CALC_DIR), - "order": 2.0, - }, - # "Number of threads": { - # "type": "integer", - # #"label": "Number of cores", - # "minimum": 1, - # "default": 1, - # "order": 3.0 - # }, - # "Memory per core": { - # "type": "integer", - # #"label" "Memory per core", - # "minimum": 1, - # "default": 1, - # "suffix": " GB", - # "order": 4.0 - # }, - "help": { - "type": "text", - "label": "For help see", - "default": "https://crest-lab.github.io/crest-docs/", - "order": 9.0, - }, - "solvent": { - "type": "stringList", - "label": "Solvation", - "values": [ - "none", - "acetone", - "acetonitrile", - "aniline", - "benzaldehyde", - "benzene", - "ch2cl2", - "chcl3", - "cs2", - "dioxane", - "dmf", - "dmso", - "ether", - "ethylacetate", - "furane", - "hexandecane", - "hexane", - "methanol", - "nitromethane", - "octanol", - "woctanol", - "phenol", - "toluene", - "thf", - "water", - ], - "default": 0, - "order": 6.0, - }, "ewin": { "type": "integer", "label": "Keep all conformers within", @@ -114,15 +47,18 @@ "default": False, "order": 8.0, }, + "help": { + "type": "text", + "label": "For help see", + "default": "https://crest-lab.github.io/crest-docs/", + "order": 9.0, + }, }, } # Display energy in kcal if user has insisted on it if easyxtb.config["energy_units"] == "kcal/mol": options["userOptions"]["ewin"]["default"] = 6 options["userOptions"]["ewin"]["suffix"] = " kcal/mol" - # Make solvation default if found in user config - if easyxtb.config["solvent"] is not None: - options["userOptions"]["solvent"]["default"] = easyxtb.config["solvent"] print(json.dumps(options)) if args.display_name: print("Conformers…") @@ -135,30 +71,17 @@ # Extract the coords geom = easyxtb.Geometry.from_cjson(avo_input["cjson"]) - # If provided crest path different to that stored, use it and save it - if Path(avo_input["crest_bin"]) != easyxtb.CREST_BIN: - crest_bin = Path(avo_input["crest_bin"]) - easyxtb.config["crest_bin"] = str(crest_bin) - with open(easyxtb.config_file, "w", encoding="utf-8") as config_path: - json.dump(easyxtb.config, config_path) - # crest takes energies in kcal so convert if provided in kJ (default) if easyxtb.config["energy_units"] == "kJ/mol": ewin_kcal = avo_input["ewin"] / 4.184 else: ewin_kcal = avo_input["ewin"] - - # Convert "none" string to Python None - if avo_input["solvent"] == "none": - solvation = None - else: - solvation = avo_input["solvent"] # Run calculation; returns set of conformers as well as Calculation object conformers, calc = easyxtb.calculate.conformers( geom, - solvation=solvation, - method=2, + solvation=easyxtb.config["solvent"], + method=easyxtb.config["method"], ewin=ewin_kcal, hess=avo_input["hess"], return_calc=True, @@ -190,13 +113,6 @@ with open(easyxtb.TEMP_DIR / "result.cjson", "w", encoding="utf-8") as save_file: json.dump(output["cjson"], save_file, indent=2) - # If user specified a save location, copy calculation directory to there - if not ( - avo_input["save_dir"] in ["", None] - or Path(avo_input["save_dir"]) == easyxtb.TEMP_DIR - ): - copytree(easyxtb.TEMP_DIR, Path(avo_input["save_dir"]), dirs_exist_ok=True) - # Pass back to Avogadro print(json.dumps(output)) logger.debug(f"The following dictionary was passed back to Avogadro: {output}") diff --git a/avo_xtb/deprotonate.py b/avo_xtb/deprotonate.py index 06d93bb..7f172ba 100644 --- a/avo_xtb/deprotonate.py +++ b/avo_xtb/deprotonate.py @@ -45,7 +45,7 @@ tautomers, calc = easyxtb.calculate.deprotonate( geom, solvation=easyxtb.config["solvent"], - method=2, + method=easyxtb.config["method"], return_calc=True, ) diff --git a/avo_xtb/protonate.py b/avo_xtb/protonate.py index 5f33021..3f02aab 100644 --- a/avo_xtb/protonate.py +++ b/avo_xtb/protonate.py @@ -5,6 +5,7 @@ import json import logging import sys +from copy import deepcopy from support import easyxtb @@ -19,7 +20,7 @@ def cleanup_after_taut(cjson: dict) -> dict: Essentially gives an empty cjson, with only the total charge and spin retained. """ - output = easyxtb.convert.empty_cjson.copy() + output = deepcopy(easyxtb.convert.empty_cjson) output["properties"]["totalCharge"] = cjson["properties"]["totalCharge"] output["properties"]["totalSpinMultiplicity"] = cjson["properties"]["totalSpinMultiplicity"] diff --git a/avo_xtb/run_xtb.py b/avo_xtb/run_xtb.py index ddc39b2..a64dcb5 100644 --- a/avo_xtb/run_xtb.py +++ b/avo_xtb/run_xtb.py @@ -7,6 +7,7 @@ import json import logging import sys +from copy import deepcopy from pathlib import Path from shutil import copytree @@ -120,7 +121,7 @@ # Start by passing back an empty cjson, then add changes output = { "moleculeFormat": "cjson", - "cjson": easyxtb.convert.empty_cjson.copy(), + "cjson": deepcopy(easyxtb.convert.empty_cjson), } # TODO Catch errors in xtb execution diff --git a/avo_xtb/solvate.py b/avo_xtb/solvate.py new file mode 100644 index 0000000..d2273d6 --- /dev/null +++ b/avo_xtb/solvate.py @@ -0,0 +1,181 @@ +# SPDX-FileCopyrightText: 2024 Matthew J. Milner +# SPDX-License-Identifier: BSD-3-Clause + +import argparse +import json +import logging +import sys +from copy import deepcopy +from pathlib import Path + +from support import easyxtb +from about import XTB_VERSION, CREST_VERSION + + +logger = logging.getLogger(__name__) + + +def split_cjson_by_layer(cjson: dict) -> list[dict]: + """Separate a CJSON into multiple CJSONs according to which layer atoms are in. + + For now drops all data except atom and bond information. + + All layers are assumed to contain neutral singlets, except layer 0, which is given + the total charge and multiplicity of the original CJSON. + """ + output = [] + for layer_index in range(max(cjson["atoms"]["layer"]) + 1): + layer = deepcopy(easyxtb.convert.empty_cjson) + atoms_in_layer = [] + for atom_index, atom_layer in enumerate(cjson["atoms"]["layer"]): + if atom_layer == layer_index: + atoms_in_layer.append(atom_index) + for atom_index in atoms_in_layer: + layer["atoms"]["coords"]["3d"].extend( + cjson["atoms"]["coords"]["3d"][3*atom_index:3*(atom_index+1)] + ) + layer["atoms"]["elements"]["number"].append( + cjson["atoms"]["elements"]["number"][atom_index] + ) + for bond_index, bond_order in enumerate(cjson["bonds"]["order"]): + bond_members = ( + cjson["bonds"]["connections"]["index"][2*bond_index:2*(bond_index+1)] + ) + if all([atom in atoms_in_layer for atom in bond_members]): + layer["bonds"]["connections"]["index"].extend(bond_members) + layer["bonds"]["order"].append(bond_order) + if layer_index == 0: + layer["properties"]["totalCharge"] = cjson["properties"]["totalCharge"] + layer["properties"]["totalSpinMultiplicity"] = cjson["properties"]["totalSpinMultiplicity"] + output.append(layer) + return output + + +def solvate(avo_input: dict) -> dict: + # Extract the cjson + full_cjson = avo_input["cjson"] + # Sort atoms based on layer + layers = split_cjson_by_layer(full_cjson) + + # Adjust for difference in indexing between what the user's choice was based on + # (Avogadro GUI uses 1-indexing) and what we receive (CJSON uses 0-indexing) + solute_layer = avo_input["solute_layer"] - 1 + solvent_layer = avo_input["solvent_layer"] - 1 + + solute_cjson = layers[solute_layer] + solvent_cjson = layers[solvent_layer] + solute_geom = easyxtb.Geometry.from_cjson(solute_cjson) + solvent_geom = easyxtb.Geometry.from_cjson(solvent_cjson) + + # Run calculation; returns new Geometry + output_geom = easyxtb.calculate.solvate( + solute_geom, + solvent_geom, + nsolv=avo_input["nsolv"], + method=2, + ) + + # Add solute bonding information from input CJSON + output_cjson = output_geom.to_cjson() + output_cjson["bonds"] = deepcopy(solute_cjson["bonds"]) + + # Format everything appropriately for Avogadro + output = { + "moleculeFormat": "cjson", + "cjson": output_cjson, + } + + # Save result + with open(easyxtb.TEMP_DIR / "result.cjson", "w", encoding="utf-8") as save_file: + json.dump(output["cjson"], save_file, indent=2) + + return output + + +if __name__ == "__main__": + parser = argparse.ArgumentParser() + parser.add_argument("--debug", action="store_true") + parser.add_argument("--print-options", action="store_true") + parser.add_argument("--run-command", action="store_true") + parser.add_argument("--display-name", action="store_true") + parser.add_argument("--lang", nargs="?", default="en") + parser.add_argument("--menu-path", action="store_true") + parser.add_argument("--test", action="store_true") + args = parser.parse_args() + + # Disable if xtb or crest missing + if easyxtb.XTB_BIN is None or easyxtb.CREST_BIN is None: + quit() + + if args.print_options: + options = { + "userOptions": { + "solute_layer": { + "type": "integer", + "label": "Layer containing the solute", + "default": 1, + "minimum": 1, + "maximum": 100, + "order": 1.0, + }, + "solvent_layer": { + "type": "integer", + "label": "Layer containing a solvent molecule", + "default": 2, + "minimum": 1, + "maximum": 100, + "order": 2.0, + }, + "nsolv": { + "type": "integer", + "label": "Number of solvent molecules to add", + "default": 4, + "minimum": 1, + "maximum": 1000, + "order": 3.0, + }, + }, + } + # If we think it won't work with the versions of xtb and CREST being used, show + # a warning + xtb_version_tuple = tuple([int(n) for n in XTB_VERSION.split(".")]) + crest_version_tuple = tuple([int(n) for n in CREST_VERSION.split(".")]) + if crest_version_tuple >= (3,0,0) and xtb_version_tuple < (6,7,0): + options["userOptions"]["warning"] = { + "type": "text", + "label": "WARNING", + "default": "QCG calculations with version 3 of CREST require xtb 6.7.1 or newer, so this command is unlikely to work for you!", + "order": 4.0, + } + elif (3,0,0) <= crest_version_tuple <= (3,1,0): + options["userOptions"]["warning"] = { + "type": "text", + "label": "WARNING", + "default": "QCG calculations with versions 3.0.x of CREST have a bug, so this command is unlikely to work for you!", + "order": 4.0, + } + print(json.dumps(options)) + + if args.display_name: + print("Solvate…") + + if args.menu_path: + print("Extensions|Semi-empirical (xtb){760}") + + if args.run_command: + # Read input from Avogadro + avo_input = json.loads(sys.stdin.read()) + output = solvate(avo_input) + # Pass back to Avogadro + print(json.dumps(output)) + logger.debug(f"The following dictionary was passed back to Avogadro: {output}") + + if args.test: + with open( + Path(__file__).parent / "tests/files/solvate.cjson", encoding="utf-8" + ) as f: + avo_input = json.load(f) + cjson = avo_input["cjson"] + layers = split_cjson_by_layer(cjson) + print(layers[0]) + print(layers[1]) diff --git a/avo_xtb/tests/files/solvate.cjson b/avo_xtb/tests/files/solvate.cjson new file mode 100644 index 0000000..41d8bc2 --- /dev/null +++ b/avo_xtb/tests/files/solvate.cjson @@ -0,0 +1,269 @@ +{ + "charge": 0, + "cjson": { + "atoms": { + "coords": { + "3d": [ + -3.44253755444722, + -0.28076706927236, + 0.08202211513711, + -2.19619354155911, + 0.5377535234173, + -0.14991040912003, + -3.30959216565811, + -0.96499512043363, + 0.91550743984568, + -4.29385343018696, + 0.37029680791868, + 0.26126111819385, + -3.6393686864466, + -0.86858435111878, + -0.81356324764349, + -2.22708173547159, + 1.70867650958983, + -0.42661735130905, + -0.90939274295079, + -0.23608552578054, + -0.00534584276377, + -0.07606115424566, + 0.33873568130699, + -0.39972676900215, + -0.73355640208978, + -0.42786150599808, + 1.05196261884314, + -0.97902512498532, + -1.19297157140393, + -0.51587626847198, + -3.6038336753845215, + 3.6284279823303223, + 0.5011339783668518, + -3.1711597442626953, + 3.5739572048187256, + 1.5878840684890747, + -3.159013271331787, + 2.743358612060547, + -0.09939655661582947 + ] + }, + "elements": { + "number": [ + 6, + 6, + 1, + 1, + 1, + 8, + 6, + 1, + 1, + 1, + 8, + 1, + 1 + ] + }, + "formalCharges": [ + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0 + ], + "layer": [ + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 1, + 1, + 1 + ] + }, + "bonds": { + "connections": { + "index": [ + 0, + 1, + 0, + 2, + 0, + 3, + 0, + 4, + 1, + 5, + 1, + 6, + 6, + 7, + 6, + 8, + 6, + 9, + 10, + 11, + 10, + 12 + ] + }, + "order": [ + 1, + 1, + 1, + 1, + 2, + 1, + 1, + 1, + 1, + 1, + 1 + ] + }, + "chemicalJson": 1, + "layer": { + "enable": { + "Ball and Stick": [ + true, + true + ], + "Cartoons": [ + true, + true + ], + "Close Contacts": [ + false, + false + ], + "Crystal Lattice": [ + false, + false + ], + "Force": [ + false, + false + ], + "Labels": [ + false, + false + ], + "Licorice": [ + false, + false + ], + "Meshes": [ + false, + false + ], + "Non-Covalent": [ + false, + false + ], + "Van der Waals": [ + false, + false + ], + "Wireframe": [ + false, + false + ] + }, + "locked": [ + false, + false, + false + ], + "settings": { + "Ball and Stick": [ + "true true 0.300000 0.100000", + "true true 0.300000 0.100000" + ], + "Cartoons": [ + "false false false false false true false", + "false false false false false true false" + ] + }, + "visible": [ + true, + true, + true + ] + }, + "properties": { + "fileName": "/home/matt/Desktop/layers.cjson", + "modelView": [ + [ + 0.72907954454422, + -0.677324652671814, + 0.09824412316083908, + -0.8313446044921875 + ], + [ + -0.24637705087661743, + -0.3936602771282196, + -0.8856112957000732, + 0.5570162534713745 + ], + [ + 0.6385230422019958, + 0.6214767694473267, + -0.4538900554180145, + -11.671831130981445 + ], + [ + 0, + 0, + 0, + 1 + ] + ], + "projection": [ + [ + 1.112694263458252, + 0, + 0, + 0 + ], + [ + 0, + 2, + 0, + 0 + ], + [ + 0, + 0, + -1.2318322658538818, + -4.463664531707764 + ], + [ + 0, + 0, + -1, + 0 + ] + ], + "totalCharge": 0, + "totalSpinMultiplicity": 1 + } + }, + "nsolv": 4, + "selectedAtoms": [], + "solute_layer": 1, + "solvent_layer": 2, + "spin": 1 +} \ No newline at end of file diff --git a/easyxtb/src/easyxtb/calculate.py b/easyxtb/src/easyxtb/calculate.py index a395bbe..a392a13 100644 --- a/easyxtb/src/easyxtb/calculate.py +++ b/easyxtb/src/easyxtb/calculate.py @@ -75,9 +75,9 @@ def __init__( self, program: str | os.PathLike = "xtb", runtype: str | None = None, - runtype_args: list[str] | None = None, + runtype_args: list | None = None, options: dict | None = None, - command: list[str] | None = None, + command: list | None = None, input_geometry: Geometry | None = None, calc_dir: os.PathLike | None = None, ): @@ -90,6 +90,10 @@ def __init__( To use a flag that takes no argument, set the value to `True`. Flags with values of `False` or `None` will not be passed. + `Geometry` objects passed as items in `runtype_args` or `command` or as values + in `options` will be saved as XYZ files and the path to the file will be put in + the command in its place. + Note that any contents of `calc_dir` will be removed when the calculation begins. """ @@ -127,6 +131,79 @@ def __init__( self.calc_dir = Path(calc_dir) if calc_dir else TEMP_DIR logger.info(f"New Calculation created for runtype {self.runtype} with {self.program}") + def _build_xtb_command(self, geom_file): + # Build command line args + # "xtb" + command = [self.program_path] + # Charge and spin from the initial Geometry + if "chrg" not in self.options and self.input_geometry.charge != 0: + command.extend(["--chrg", self.input_geometry.charge]) + if "uhf" not in self.options and self.input_geometry.spin != 0: + command.extend(["--uhf", self.input_geometry.spin]) + # Any other options + for flag, value in self.options.items(): + # Add appropriate number of minuses to flags + if len(flag) == 1: + flag = "-" + flag + else: + flag = "--" + flag + if value is True: + command.append(flag) + elif value is False or value is None: + continue + elif isinstance(value, Geometry): + command.extend([flag, value]) + else: + command.extend([flag, value]) + # Which calculation to run + if self.runtype is None: + # Simple single point + pass + else: + command.extend([ + "--" + self.runtype, + *self.runtype_args, + ]) + # Path to the input geometry file, preceded by a -- to ensure that it is not + # parsed as an argument to the runtype + command.extend(["--", geom_file]) + return command + + def _build_crest_command(self, geom_file): + # Build command line args + # "crest" + command = [self.program_path] + # Path to the input geometry file + command.append(geom_file) + # Which calculation to run + if self.runtype is None: + # v3 iMTD-GC sampling + pass + else: + command.extend([ + "--" + self.runtype, + *self.runtype_args, + ]) + # Charge and spin from the initial Geometry + if "chrg" not in self.options and self.input_geometry.charge != 0: + command.extend(["--chrg", self.input_geometry.charge]) + if "uhf" not in self.options and self.input_geometry.spin != 0: + command.extend(["--uhf", self.input_geometry.spin]) + # Any other options + for flag, value in self.options.items(): + # Add appropriate number of minuses to flags + if len(flag) == 1: + flag = "-" + flag + else: + flag = "--" + flag + if value is True: + command.append(flag) + elif value is False or value is None: + continue + else: + command.extend([flag, value]) + return command + def run(self): """Run calculation with xtb or crest, storing the output, the saved output file, and the parsed energy, as well as the `subprocess` object.""" @@ -145,7 +222,7 @@ def run(self): self.calc_dir.mkdir(parents=True) # Save geometry to file - geom_file = self.calc_dir / "input.xyz" + geom_file = self.calc_dir/"input.xyz" logger.debug(f"Saving input geometry to {geom_file}") self.input_geometry.write_xyz(geom_file) # We are using proper paths for pretty much everything so it shouldn't be @@ -156,44 +233,24 @@ def run(self): logger.debug(f"Working directory changed to {Path.cwd()}") self.output_file = geom_file.with_name("output.out") + # Build command line args if self.command: # If arguments were passed by the user, use them as is command = self.command + elif self.program == "crest": + # CREST parses command line input in a slightly different way to xtb + command = self._build_crest_command(geom_file) else: - # Build command line args - # "xtb" - command = [str(self.program_path)] - # Charge and spin from the initial Geometry - if "chrg" not in self.options and self.input_geometry.charge != 0: - command.extend(["--chrg", str(self.input_geometry.charge)]) - if "uhf" not in self.options and self.input_geometry.spin != 0: - command.extend(["--uhf", str(self.input_geometry.spin)]) - # Any other options - for flag, value in self.options.items(): - # Add appropriate number of minuses to flags - if len(flag) == 1: - flag = "-" + flag - else: - flag = "--" + flag - if value is True: - command.append(flag) - elif value is False or value is None: - continue - else: - command.extend([flag, str(value)]) - # Which calculation to run - if self.runtype is None: - # Simple single point - pass - else: - command.extend([ - "--" + self.runtype, - *self.runtype_args, - "--", - ]) - # Path to the input geometry file, preceded by a -- to ensure that it is not - # parsed as an argument to the runtype - command.extend(["--", str(geom_file)]) + command = self._build_xtb_command(geom_file) + # Replace any Geometry objects with paths to files + aux_count = 0 + for i, arg in enumerate(command): + if isinstance(arg, Geometry): + aux_count += 1 + aux_file = arg.write_xyz(self.calc_dir/f"aux{aux_count}.xyz") + command[i] = aux_file + # Sanitize everything to strings + command = [x if isinstance(x, str) else str(x) for x in command] logger.debug(f"Calculation will be run with the command: {' '.join(command)}") # Run xtb or crest from command line @@ -218,7 +275,6 @@ def run(self): self.process_xtb() def process_xtb(self): - # First do generic operations that apply to many calculation types # Extract energy from output (if not found, returns None) self.energy = parse_energy(self.output) @@ -332,6 +388,16 @@ def process_crest(self): self.output_geometry = best["geometry"] self.energy = best["energy"] logger.debug(f"Energy of lowest energy tautomer: {self.energy}") + case "qcg": + # Explicit solvent shell growing + logger.debug("Growing of a solvent shell was requested, so checking for generated cluster geometry") + # Get final cluster geom + self.output_geometry = Geometry.from_file( + self.calc_dir/"grow/cluster.xyz", + charge=self.input_geometry.charge, + spin=self.input_geometry.spin, + ) + logger.debug(f'Cluster geometry read from {self.calc_dir/"grow/cluster.xyz"}') case _: pass @@ -504,7 +570,7 @@ def conformers( hess: bool = False, n_proc: int | None = None, return_calc: bool = False, -) -> list[dict]: +) -> list[dict] | tuple[list[dict], Calculation]: """Simulate a conformer ensemble and return set of conformer Geometries and energies. The returned conformers are ordered from lowest to highest energy. @@ -540,7 +606,7 @@ def protonate( method: int = 2, n_proc: int | None = None, return_calc: bool = False, -) -> list[dict]: +) -> list[dict] | tuple[list[dict], Calculation]: """Screen possible protonation sites and return set of tautomer Geometries and energies. The returned tautomers are ordered from lowest to highest energy. @@ -570,7 +636,7 @@ def deprotonate( method: int = 2, n_proc: int | None = None, return_calc: bool = False, -) -> list[dict]: +) -> list[dict] | tuple[list[dict], Calculation]: """Screen possible deprotonation sites and return set of tautomer Geometries and energies. The returned tautomers are ordered from lowest to highest energy. @@ -592,3 +658,37 @@ def deprotonate( return calc.tautomers, calc else: return calc.tautomers + + +def solvate( + solute_geometry: Geometry, + solvent_geometry: Geometry, + nsolv: int, + method: int = 2, + n_proc: int | None = None, + return_calc: bool = False, +) -> Geometry | tuple[Geometry, Calculation]: + """Grow a solvent shell around a solute for a total of `nsolv` solvent molecules. + + Note that non-zero charge and spin on the solvent `Geometry` will not be passed to + CREST. + """ + method_flag = f"gfn{method}" + calc = Calculation( + program="crest", + input_geometry=solute_geometry, + runtype="qcg", + runtype_args=[solvent_geometry], + options={ + "grow": True, + "nsolv": nsolv, + "xnam": XTB_BIN, + method_flag: True, + "T": n_proc if n_proc else config["n_proc"], + }, + ) + calc.run() + if return_calc: + return calc.output_geometry, calc + else: + return calc.output_geometry diff --git a/easyxtb/src/easyxtb/convert.py b/easyxtb/src/easyxtb/convert.py index 31b3c78..cc5c823 100644 --- a/easyxtb/src/easyxtb/convert.py +++ b/easyxtb/src/easyxtb/convert.py @@ -42,7 +42,6 @@ def convert_freq(freq=None, wavelength=None, wavenumber=None): "number": [], }, "formalCharges": [], - "layer": [], }, "bonds": { "connections": { diff --git a/plugin.json b/plugin.json index d30764f..b4d1f1e 100644 --- a/plugin.json +++ b/plugin.json @@ -20,6 +20,7 @@ { "name": "conformers", "command": "avo_xtb/conformers.py" }, { "name": "protonate", "command": "avo_xtb/protonate.py" }, { "name": "deprotonate", "command": "avo_xtb/deprotonate.py" }, + { "name": "solvate", "command": "avo_xtb/solvate.py" }, { "name": "crest-docs", "command": "avo_xtb/crest_docs.py" } ] }