diff --git a/CHANGELOG.md b/CHANGELOG.md index e5ba49e..f46ded0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,10 +7,15 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] ### Changed - vdW radii scaling parameter can now be adjusted via `mindlessgen.toml` or CLI +- better type hints for `Callables` ### Fixed - Unit conversion for (currenly unused) vdW radii from the original Fortran project -- minor print output issues (no new line breaks...) +- minor print output issues (no new line breaks, more consistent verbosity differentiation, ...) +- bug in `postprocess_mol` which led to an unassigned return variable in the single-point case + +### Added +- Support for the novel "g-xTB" method (working title: GP3-xTB) ## [0.4.0] - 2024-09-19 ### Changed diff --git a/src/mindlessgen/generator/main.py b/src/mindlessgen/generator/main.py index e8d4c13..237d1ec 100644 --- a/src/mindlessgen/generator/main.py +++ b/src/mindlessgen/generator/main.py @@ -4,12 +4,13 @@ from __future__ import annotations +from collections.abc import Callable from pathlib import Path import multiprocessing as mp import warnings from ..molecules import generate_random_molecule, Molecule -from ..qm import XTB, get_xtb_path, QMMethod, ORCA, get_orca_path +from ..qm import XTB, get_xtb_path, QMMethod, ORCA, get_orca_path, GP3, get_gp3_path from ..molecules import iterative_optimization, postprocess_mol from ..prog import ConfigManager @@ -41,12 +42,15 @@ def generator(config: ConfigManager) -> tuple[list[Molecule] | None, int]: # Import and set up required engines refine_engine: QMMethod = setup_engines( - config.refine.engine, config, get_xtb_path, get_orca_path + config.refine.engine, + config, + get_xtb_path, + get_orca_path, # GP3 cannot be used anyway ) if config.general.postprocess: postprocess_engine: QMMethod | None = setup_engines( - config.postprocess.engine, config, get_xtb_path, get_orca_path + config.postprocess.engine, config, get_xtb_path, get_orca_path, get_gp3_path ) else: postprocess_engine = None @@ -273,8 +277,9 @@ def header(version: str) -> str: def setup_engines( engine_type: str, cfg: ConfigManager, - xtb_path_func, - orca_path_func, + xtb_path_func: Callable, + orca_path_func: Callable, + gp3_path_func: Callable | None = None, ): """ Set up the required engine. @@ -295,5 +300,12 @@ def setup_engines( except ImportError as e: raise ImportError("orca not found.") from e return ORCA(path, cfg.orca) + elif engine_type == "gp3": + if gp3_path_func is None: + raise ImportError("No callable function for determining the gp3 path.") + path = gp3_path_func() + if not path: + raise ImportError("'gp3' binary could not be found.") + return GP3(path) else: raise NotImplementedError("Engine not implemented.") diff --git a/src/mindlessgen/molecules/postprocess.py b/src/mindlessgen/molecules/postprocess.py index 543285d..17f6ae4 100644 --- a/src/mindlessgen/molecules/postprocess.py +++ b/src/mindlessgen/molecules/postprocess.py @@ -34,6 +34,7 @@ def postprocess_mol( else: try: engine.singlepoint(mol, verbosity=verbosity) + postprocmol = mol except RuntimeError as e: raise RuntimeError( "Single point calculation in postprocessing failed." diff --git a/src/mindlessgen/molecules/refinement.py b/src/mindlessgen/molecules/refinement.py index 638df9b..291c1a9 100644 --- a/src/mindlessgen/molecules/refinement.py +++ b/src/mindlessgen/molecules/refinement.py @@ -68,7 +68,7 @@ def iterative_optimization( ) break # Stop if the atom counts are the same as in the previous cycle if len(fragmols) == 1: - if verbosity > 0: + if verbosity > 1: print( f"Only one fragment detected in cycle {cycle + 1}. " + f"Stopping at cycle {cycle + 1}." diff --git a/src/mindlessgen/prog/config.py b/src/mindlessgen/prog/config.py index 6562729..491b4b7 100644 --- a/src/mindlessgen/prog/config.py +++ b/src/mindlessgen/prog/config.py @@ -497,7 +497,7 @@ def engine(self, engine: str): """ if not isinstance(engine, str): raise TypeError("Postprocess engine should be a string.") - if engine not in ["xtb", "orca"]: + if engine not in ["xtb", "orca", "gp3"]: raise ValueError("Postprocess engine can only be xtb or orca.") self._engine = engine diff --git a/src/mindlessgen/qm/__init__.py b/src/mindlessgen/qm/__init__.py index 18cc01a..a5e8e3a 100644 --- a/src/mindlessgen/qm/__init__.py +++ b/src/mindlessgen/qm/__init__.py @@ -5,5 +5,14 @@ from .base import QMMethod from .xtb import XTB, get_xtb_path from .orca import ORCA, get_orca_path +from .gp3 import GP3, get_gp3_path -__all__ = ["XTB", "get_xtb_path", "QMMethod", "ORCA", "get_orca_path"] +__all__ = [ + "XTB", + "get_xtb_path", + "QMMethod", + "ORCA", + "get_orca_path", + "GP3", + "get_gp3_path", +] diff --git a/src/mindlessgen/qm/gp3.py b/src/mindlessgen/qm/gp3.py new file mode 100644 index 0000000..0ca4d40 --- /dev/null +++ b/src/mindlessgen/qm/gp3.py @@ -0,0 +1,180 @@ +""" +This module contains all interactions with the GP3-xTB binary +for next-gen tight-binding calculations. +""" + +import subprocess as sp +import shutil +from pathlib import Path +from tempfile import TemporaryDirectory + +from ..molecules import Molecule +from .base import QMMethod + + +class GP3(QMMethod): + """ + This class handles all interaction with the GP3 external dependency. + """ + + def __init__(self, path: str | Path) -> None: + """ + Initialize the GP3 class. + """ + if isinstance(path, str): + self.path: Path = Path(path).resolve() + elif isinstance(path, Path): + self.path = path + else: + raise TypeError("gp3_path should be a string or a Path object.") + + def singlepoint(self, molecule: Molecule, verbosity: int = 1) -> str: + """ + Perform a single-point calculation using GP3-xTB. + """ + + # Create a unique temporary directory using TemporaryDirectory context manager + with TemporaryDirectory(prefix="gp3_") as temp_dir: + temp_path = Path(temp_dir).resolve() + # write the molecule to a temporary file + molecule.write_xyz_to_file(str(temp_path / "molecule.xyz")) + + # run gp3 + arguments = [ + "-c", + "molecule.xyz", + ] + # dump molecule.charge and molecule.uhf to a file + if molecule.charge != 0: + with open(temp_path / ".CHRG", "w", encoding="utf-8") as f: + f.write(str(molecule.charge)) + if molecule.uhf != 0: + with open(temp_path / ".UHF", "w", encoding="utf-8") as f: + f.write(str(molecule.uhf)) + + if verbosity > 2: + print(f"Running command: gp3 {' '.join(arguments)}") + + gp3_log_out, gp3_log_err, return_code = self._run( + temp_path=temp_path, arguments=arguments + ) + if verbosity > 2: + print(gp3_log_out) + if return_code != 0: + raise RuntimeError( + f"GP3-xTB failed with return code {return_code}:\n{gp3_log_err}" + ) + + return gp3_log_out + + def check_gap( + self, molecule: Molecule, threshold: float = 0.5, verbosity: int = 1 + ) -> bool: + """ + Check if the HL gap is larger than a given threshold. + + Arguments: + molecule (Molecule): Molecule to check + threshold (float): Threshold for the gap + + Returns: + bool: True if the gap is larger than the threshold, False otherwise + """ + + # Perform a single point calculation + try: + gp3_out = self.singlepoint(molecule) + except RuntimeError as e: + raise RuntimeError("Single point calculation failed.") from e + + # Parse the output to get the gap + hlgap = None + for line in gp3_out.split("\n"): + if "gap (eV)" in line and "dE" not in line: + # check if "alpha->alpha" is present in the same line + # then, the line looks as follows: + # gap (eV) alpha->alpha : 7.02263204 + if "alpha->alpha" in line: + hlgap = float(line.split()[4]) + break + # otherwise, the line looks as follows: + # gap (eV) : 9.06694119 + hlgap = float(line.split()[3]) + break + if hlgap is None: + raise ValueError("GP3-xTB gap not determined.") + + if verbosity > 1: + print(f"GP3-xTB HOMO-LUMO gap: {hlgap:5f}") + + return hlgap > threshold + + def optimize( + self, molecule: Molecule, max_cycles: int | None = None, verbosity: int = 1 + ) -> Molecule: + """ + Optimize a molecule using GP3-xTB. + """ + raise NotImplementedError("Optimization is not yet implemented for GP3-xTB.") + + def _run(self, temp_path: Path, arguments: list[str]) -> tuple[str, str, int]: + """ + Run GP3-xTB with the given arguments. + + Arguments: + arguments (list[str]): The arguments to pass to GP3-xTB. + + Returns: + tuple[str, str, int]: The output of the GP3-xTB calculation (stdout and stderr) + and the return code + """ + try: + gp3_out = sp.run( + [str(self.path)] + arguments, + cwd=temp_path, + capture_output=True, + check=True, + ) + # get the output of the GP3-xTB calculation (of both stdout and stderr) + gp3_log_out = gp3_out.stdout.decode("utf8") + gp3_log_err = gp3_out.stderr.decode("utf8") + if ( + "no SCF convergence" in gp3_log_out + or "nuclear repulsion" not in gp3_log_out + ): + raise sp.CalledProcessError( + 1, + str(self.path), + gp3_log_out.encode("utf8"), + gp3_log_err.encode("utf8"), + ) + return gp3_log_out, gp3_log_err, 0 + except sp.CalledProcessError as e: + gp3_log_out = e.stdout.decode("utf8") + gp3_log_err = e.stderr.decode("utf8") + return gp3_log_out, gp3_log_err, e.returncode + + +# TODO: 1. Convert this to a @staticmethod of Class GP3 +# 2. Rename to `get_method` or similar to enable an abstract interface +# 3. Add the renamed method to the ABC `QMMethod` +# 4. In `main.py`: Remove the passing of the path finder functions as arguments +# and remove the boiler plate code to make it more general. +def get_gp3_path(binary_name: str | Path | None = None) -> Path: + """ + Get the path to the GP3 binary based on different possible names + that are searched for in the PATH. + """ + default_gp3_names: list[str | Path] = ["gp3", "gp3_dev"] + # put binary name at the beginning of the lixt to prioritize it + if binary_name is not None: + binary_names = [binary_name] + default_gp3_names + else: + binary_names = default_gp3_names + # Get gp3 path from 'which gp3' command + for binpath in binary_names: + which_gp3 = shutil.which(binpath) + if which_gp3: + gp3_path = Path(which_gp3).resolve() + return gp3_path + raise ImportError("'gp3' binary could not be found.") diff --git a/src/mindlessgen/qm/xtb.py b/src/mindlessgen/qm/xtb.py index 163b936..60d3afd 100644 --- a/src/mindlessgen/qm/xtb.py +++ b/src/mindlessgen/qm/xtb.py @@ -76,7 +76,7 @@ def optimize( def singlepoint(self, molecule: Molecule, verbosity: int = 1) -> str: """ - Optimize a molecule using xtb. + Perform a single-point calculation using xtb. """ # Create a unique temporary directory using TemporaryDirectory context manager @@ -97,7 +97,7 @@ def singlepoint(self, molecule: Molecule, verbosity: int = 1) -> str: arguments += ["--uhf", str(molecule.uhf)] if verbosity > 2: - print(f"Running command: {' '.join(arguments)}") + print(f"Running command: xtb {' '.join(arguments)}") xtb_log_out, xtb_log_err, return_code = self._run( temp_path=temp_path, arguments=arguments @@ -141,7 +141,7 @@ def check_gap( if hlgap is None: raise ValueError("HOMO-LUMO gap not determined.") if verbosity > 1: - print(f"HOMO-LUMO gap: {hlgap:5f}") + print(f"xTB HOMO-LUMO gap: {hlgap:5f}") if hlgap > threshold: return True