Skip to content

Commit

Permalink
Add support for the new "g-xTB" method (working title: GP3-xTB) and m…
Browse files Browse the repository at this point in the history
…inor bug fixes (grimme-lab#49)

* placeholder for GP3-xTB implementation

Signed-off-by: Marcel Müller <marcel.mueller@thch.uni-bonn.de>

* implementation of GP3 as a QMMethod

Signed-off-by: Marcel Mueller <marcel.mueller@thch.uni-bonn.de>

* fix for picking wrong gap(ev) line and better printout

Signed-off-by: Marcel Mueller <marcel.mueller@thch.uni-bonn.de>

* update CHANGELOG.md

Signed-off-by: Marcel Müller <marcel.mueller@thch.uni-bonn.de>

---------

Signed-off-by: Marcel Müller <marcel.mueller@thch.uni-bonn.de>
Signed-off-by: Marcel Mueller <marcel.mueller@thch.uni-bonn.de>
marcelmbn authored Sep 27, 2024
1 parent 376fe6b commit 5a38779
Showing 8 changed files with 219 additions and 12 deletions.
7 changes: 6 additions & 1 deletion CHANGELOG.md
Original file line number Diff line number Diff line change
@@ -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
22 changes: 17 additions & 5 deletions src/mindlessgen/generator/main.py
Original file line number Diff line number Diff line change
@@ -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.")
1 change: 1 addition & 0 deletions src/mindlessgen/molecules/postprocess.py
Original file line number Diff line number Diff line change
@@ -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."
2 changes: 1 addition & 1 deletion src/mindlessgen/molecules/refinement.py
Original file line number Diff line number Diff line change
@@ -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}."
2 changes: 1 addition & 1 deletion src/mindlessgen/prog/config.py
Original file line number Diff line number Diff line change
@@ -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

11 changes: 10 additions & 1 deletion src/mindlessgen/qm/__init__.py
Original file line number Diff line number Diff line change
@@ -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",
]
180 changes: 180 additions & 0 deletions src/mindlessgen/qm/gp3.py
Original file line number Diff line number Diff line change
@@ -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.")
6 changes: 3 additions & 3 deletions src/mindlessgen/qm/xtb.py
Original file line number Diff line number Diff line change
@@ -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

0 comments on commit 5a38779

Please sign in to comment.