xtb

xtb: the semiempirical GFN-xTB, GFN-FF, and g-xTB program, run on the command line.

The codec writes structure.xyz and runs::

xtb structure.xyz --<runtype> --chrg <charge> --uhf <unpaired> <method flags> <keywords>

model.method selects the Hamiltonian: gfn0, gfn1, gfn2, gfnff, or gxtb. g-xTB needs the xtb build distributed with g-xTB (https://github.com/grimme-lab/g-xtb); point qccompute at it with executables={"xtb": ...} or QCCOMPUTE_XTB_EXECUTABLE. GFN0 reads its parameter file from the directory named by XTBPATH.

Keywords become command-line flags in xtb's own spelling: {"alpb": "water"} gives --alpb water, {"etemp": 1000} gives --etemp 1000, a list gives several values ({"alpb": ["water", "bar1M"]}), True gives a bare flag, and False or None leaves the flag out. For optimizations, {"opt": "tight"} sets the optimization level. Detailed input ($constrain, $fix, $wall, ...) goes in an xcontrol file supplied through files and named with {"input": "xcontrol"}; such files are passed through verbatim, so atom indices in them are 1-based as xtb expects.

Hessians also fill frequencies and normal modes from g98.out, where xtb writes mass-weighted vectors rounded to two decimals; the codec converts them to unit-length Cartesian displacements, as qcdata defines them.

One keyword is not an xtb flag: external_potentials, point charges for QM/MM electrostatic embedding, in the same form the Psi4 codec takes: [[charge, [x, y, z]], ...] in bohr. An optional third entry, [charge, [x, y, z], atomic_number], names the element whose chemical hardness xtb uses to damp that charge's interaction at short range. The codec writes them to pcharge and adds an $embedding group to the xcontrol file (appending to yours if you supply one); gradients then also fill point_charge_gradient from xtb's pcgrad. Only GFN1- and GFN2-xTB apply point charges: GFN0, GFN-FF, and g-xTB ignore them without a warning, so the codec rejects them there.

Examples

Each script runs as written. Copy it, delete the options you don't need, and run it with python script.py. They live in the repository's examples/ folder.

Single point: energy, gradient, Hessian

examples/xtb/single_point.py

"""xtb single point: energy, gradient, or Hessian with GFN-xTB, GFN-FF, or g-xTB.

Requires: ``pip install qccompute`` (it installs qccodec) and xtb on PATH, e.g.
``conda install -c conda-forge xtb``. For g-xTB (``method="gxtb"``) use the xtb build
from https://github.com/grimme-lab/g-xtb instead; it also runs every other method.
GFN0-xTB reads its parameter file from the directory named by ``XTBPATH``.

How the codec talks to xtb (https://xtb-docs.readthedocs.io/):

- ``model.method`` picks the Hamiltonian: "gfn0", "gfn1", "gfn2", "gfnff", or "gxtb".
  There is no basis set.
- Charge and multiplicity come from ``structure`` (passed as ``--chrg`` and ``--uhf``).
- ``keywords`` become xtb command-line flags (``xtb --help``): ``{"alpb": "water"}``
  gives ``--alpb water``, ``True`` gives a bare flag, and ``False``/``None`` drop it.
  Flags the codec sets itself (``chrg``, ``uhf``, ``gfn``, ``sp``, ``hess``, ...)
  are rejected.
- xtb only warns about options it does not know and runs without them; the codec
  turns that warning into a failure, so a typo never silently changes the method.
- Detailed settings (``$scc``, ``$wall``, ``$constrain``, ...) go in an xcontrol file
  (``man 7 xcontrol``); see the commented example below. Atom indices in it are
  1-based, as xtb expects.
- ``external_potentials`` (point charges, as for Psi4) is the one keyword that is not
  an xtb flag: the codec writes it to ``pcharge`` and adds an ``$embedding`` group.
- g-xTB has no ALPB/GBSA parameters (xtb stops with "No ALPB/GBSA parameters found")
  and xtb 6.7.1 silently ignores ``cpcmx`` with g-xTB, so run g-xTB in the gas phase.
- xtb's Hessian is numerical (finite differences of analytic gradients) and comes
  with harmonic frequencies and normal modes. Do not pass ``cma``: it moves the
  molecule to its principal axes, so results would not match ``structure``.
"""

from qccompute import compute
from qcdata import CalcType, Model, ProgramInput, Structure

calctype = CalcType.gradient  # or CalcType.energy, CalcType.hessian

# Coordinates are in bohr (Structure.from_xyz / Structure.open read angstrom files).
structure = Structure(
    symbols=["O", "H", "H"],
    geometry=[[0.0, 0.0, 0.0], [0.524, 1.687, 0.48], [1.146, -0.45, -1.354]],
    charge=0,
    multiplicity=1,  # open shells: xtb gets multiplicity - 1 unpaired electrons
)

prog_input = ProgramInput(
    program="xtb",
    calctype=calctype,
    structure=structure,
    model=Model(method="gfn2"),  # or "gfn1", "gfn0", "gfnff", "gxtb" (g-xTB build)
    keywords={
        # Any xtb command-line flag, without the leading dashes:
        "parallel": 1,  # threads
        "acc": 1.0,  # SCC accuracy; lower is tighter (0.0001 to 1000)
        "iterations": 250,  # maximum SCC iterations
        # "etemp": 300,  # electronic temperature (K); raise it for difficult SCCs
        #
        # Implicit solvent; pick one model (GFN1/GFN2/GFN-FF only, not g-xTB):
        "alpb": "water",  # ALPB; ["water", "bar1M"] also sets the reference state
        # "gbsa": "h2o",  # GBSA (its own solvent names)
        # "cpcmx": "water",  # CPCM-X
        #
        # Extra output, found in the logs or with collect_files=True:
        # "molden": True,  # orbitals in molden.input
        # "wbo": True,  # Wiberg bond orders in wbo
        # "json": True,  # xtbout.json
        #
        # Detailed settings from an xcontrol file (pass it in `files` below):
        # "input": "xcontrol",
        #
        # Point charges for QM/MM electrostatic embedding (GFN1/GFN2 only; xtb
        # ignores them for GFN0, GFN-FF, and g-xTB, so the codec refuses those):
        # [charge, [x, y, z]] in bohr, optionally followed by the atomic number
        # whose hardness damps the charge at short range. The gradient on the
        # charges comes back as point_charge_gradient.
        # "external_potentials": [[-0.834, [0.0, 0.0, 5.5], 8], [0.417, [1.4, 0.0, 6.4], 1]],
    },
    # files={"xcontrol": "$scc\n   temp=500\n$end\n"},
)

output = compute(
    prog_input,
    # timeout=600,  # seconds before xtb is killed
    # collect_files=True,  # keep every output file in output.results.files
    # executables={"xtb": "/path/to/g-xtb/bin/xtb"},  # or QCCOMPUTE_XTB_EXECUTABLE
)

results = output.results
print("Energy (Hartree):", results.energy)
if results.gradient is not None:
    print("Gradient (Hartree/bohr):\n", results.gradient)
if results.hessian is not None:
    print("Frequencies (cm^-1):", results.freqs_wavenumber)
    print("Normal modes (modes, atoms, xyz):", results.normal_modes_cartesian.shape)
print("Partial charges:", results.atomic_charges)
if results.point_charge_gradient is not None:
    print("Gradient on point charges:\n", results.point_charge_gradient)
print("xtb version:", results.provenance.program_version)

Optimization

examples/xtb/optimization.py

"""xtb geometry optimization with its built-in optimizer, ANCopt.

Requires: ``pip install qccompute`` (it installs qccodec) and xtb on PATH, e.g.
``conda install -c conda-forge xtb`` (or the g-xTB build of xtb for ``method="gxtb"``,
https://github.com/grimme-lab/g-xtb).

Notes for xtb optimizations (https://xtb-docs.readthedocs.io/en/latest/optimization.html):

- ``{"opt": "tight"}`` sets the convergence level: "crude", "sloppy", "loose",
  "normal" (default), "tight", "verytight", or "extreme". ``cycles`` caps the steps.
- xtb exits normally when it runs out of cycles; the codec reports that as a failure
  ("failed to converge geometry optimization").
- Each step in ``results.frames`` has its structure and energy but no gradient:
  xtb's log records only those.
- Constraints go in an xcontrol file (``$fix``, ``$constrain``, ``$wall``; see
  https://xtb-docs.readthedocs.io/en/latest/xcontrol.html). Its atom indices are
  1-based, unlike the 0-based indices used elsewhere in qcdata.
- Other keywords become xtb command-line flags, as in single_point.py: solvent,
  accuracy, electronic temperature, threads.
- For transition states, or an optimizer that records gradients at every step, drive
  xtb gradients with geomeTRIC (examples/geometric/).
"""

from qccompute import compute
from qcdata import CalcType, Model, ProgramInput, Structure

# Coordinates are in bohr.
structure = Structure(
    symbols=["O", "H", "H"],
    geometry=[[0.0, 0.0, 0.0], [0.524, 1.687, 0.48], [1.146, -0.45, -1.354]],
    charge=0,
    multiplicity=1,
)

prog_input = ProgramInput(
    program="xtb",
    calctype=CalcType.optimization,
    structure=structure,
    model=Model(method="gfn2"),  # or "gfn1", "gfnff", "gxtb" (g-xTB build)
    keywords={
        # ANCopt settings:
        "opt": "tight",  # convergence level
        "cycles": 100,  # maximum optimization steps
        #
        "parallel": 1,  # threads
        "alpb": "water",  # implicit solvent; or "gbsa": "h2o"
        # "etemp": 300,  # electronic temperature (K)
        #
        # Constraints from an xcontrol file (pass it in `files` below):
        # "input": "xcontrol",
    },
    # Keep the first O-H bond at its starting length and freeze the oxygen atom.
    # xcontrol atom indices are 1-based: atom 1 is structure.symbols[0].
    # files={
    #     "xcontrol": (
    #         "$constrain\n"
    #         "   force constant=1.0\n"
    #         "   distance: 1, 2, auto\n"
    #         "$fix\n"
    #         "   atoms: 1\n"
    #         "$end\n"
    #     )
    # },
)

output = compute(
    prog_input,
    # timeout=600,  # seconds before xtb is killed
    # collect_files=True,  # keep every output file (xtbopt.xyz, ...) in results.files
    # executables={"xtb": "/path/to/xtb"},  # or QCCOMPUTE_XTB_EXECUTABLE
)

results = output.results
print("Steps:", len(results.frames))
print("Energies (Hartree):", results.energies)
print("Final structure (bohr):\n", results.final_structure.geometry)
# Save the trajectory, e.g. to view in a molecule viewer:
# results.save("optimization.json"); results.final_structure.save("final.xyz")

Many calculations at once

examples/xtb/many_calculations.py

"""xtb for many structures at once: calculations run side by side.

Requires: what examples/xtb/single_point.py requires, and qccompute 0.21.0 or newer.

How lists work (the same for every program that runs as its own process: ORCA,
Psi4, CREST, DFTB+, TeraChem, ...):

- Pass ``compute`` a list of inputs to get a list of outputs in the same order. Each
  calculation runs in its own process and directory, so several run at once in
  threads; no ``if __name__ == "__main__":`` guard is needed.
- ``max_workers`` is a budget of CPU threads, by default the cores this process may
  use. Each calculation counts the threads it asks for: here xtb's ``parallel``
  (``-P``), for ORCA ``pal``, for Psi4 ``-n``, for CREST ``threads``. With 16 cores
  and ``parallel`` 2, eight xtb calculations run at once. GPU programs (TeraChem) get
  one GPU each through ``CUDA_VISIBLE_DEVICES``.
- The inputs may mix programs, methods, and calculation types.
- Every calculation runs first. Then, if any failed, ``ComputeErrors`` (an
  ``ExceptionGroup``) is raised; each error has the failed input's position as
  ``input_index``, and ``.outputs`` holds every output. With ``raise_exc=False`` the
  outputs are returned as they are; check ``output.success``.
- ``print_logs=True`` prints each line prefixed with its input's index, such as
  ``[3]``. An ``update_func`` receives ``(index, all_logs, new_logs)`` from the thread
  running that calculation.
"""

import math
import time

from qccompute import compute
from qccompute.exceptions import ComputeErrors
from qcdata import CalcType, Model, ProgramInput, Structure


def water(r: float) -> Structure:
    """Water with one O-H bond stretched to r bohr, the other at 1.81 bohr."""
    angle = math.radians(104.5)
    return Structure(
        symbols=["O", "H", "H"],
        geometry=[
            [0.0, 0.0, 0.0],
            [r, 0.0, 0.0],
            [1.81 * math.cos(angle), 1.81 * math.sin(angle), 0.0],
        ],
    )


distances = [1.4 + 0.1 * i for i in range(20)]  # bohr
inputs = [
    ProgramInput(
        program="xtb",
        calctype=CalcType.gradient,
        structure=water(r),
        model=Model(method="gfn2"),
        keywords={"parallel": 2},  # threads for each xtb, counted against max_workers
    )
    for r in distances
]

start = time.perf_counter()
try:
    outputs = compute(
        inputs,
        # max_workers=8,  # CPU threads all calculations may use at once
        # print_logs=True,  # stream every program's output, prefixed with [index]
        # timeout=600,  # seconds each calculation may run
    )
except ComputeErrors as errors:
    for error in errors.exceptions:
        print(f"Input {error.input_index} failed: {error}")
    outputs = errors.outputs  # every output, failed ones with success=False
parallel = time.perf_counter() - start

start = time.perf_counter()
compute(inputs, max_workers=2)  # one at a time, for comparison
serial = time.perf_counter() - start

print("r (bohr)  energy (Hartree)")
for r, output in zip(distances, outputs, strict=True):
    print(f"{r:8.2f}  {output.results.energy:16.8f}")
print(f"Side by side: {parallel:.2f} s; one at a time: {serial:.2f} s")