geomeTRIC

geomeTRIC: optimizations and transition-state searches driving another program.

geomeTRIC runs in Python and asks for gradients through a callback: every step is a full gradient calculation of the program named in subprograms, run by qccompute like any other calculation and kept in the trajectory. Any codec with gradients works as the child, including ones registered with qccodec.register.

Inputs map to geomeTRIC as follows:

• subprograms must hold one ProgramSpec with calctype="gradient": the child program, its model, and its own keywords (passed to that program unchanged).
• keywords are geomeTRIC's optimizer parameters (geometric.params.OptParams): maxiter, convergence_set, convergence_energy, trust, hessian, enforce, ... Two are handled by the codec: coordsys (tric by default, or cart, prim, dlc, hdlc) and constraints (geomeTRIC's constraint dictionary, e.g. {"freeze": [{"type": "xyz", "indices": [0]}]}, with set distances in bohr; scans are not supported).
• cmdline_args are not accepted.

Calculation types: optimization and transition_state (transition=True; geomeTRIC computes a starting Hessian from finite differences of child gradients).

The trajectory holds one frame per optimizer step. Gradients at the displaced geometries of a finite-difference Hessian are not steps, so they are left out, unless one fails. Each Hessian geomeTRIC computes is attached to the frame at its geometry (frame.results.hessian): the first frame for hessian="first", the default for transition states, and also the last frame for "first+last" or "last". Atom indices in constraints are 0-based, as in geomeTRIC's Python interface.

With qccompute, propagate_wfn=True starts each child calculation from the previous step's orbitals when the child codec supports it (TeraChem). The optimization declares the threads and GPUs its child does, so lists of optimizations share the machine correctly.

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.

Optimization

examples/geometric/optimization.py

"""geomeTRIC geometry optimization driving another program's gradients.

Requires: ``pip install qccompute qccodec[geometric]`` (geomeTRIC itself) and the
gradient program on PATH. This example uses CREST for fast GFN2-xTB gradients:
https://github.com/crest-lab/crest/releases (official static build).

How it works (https://geometric.readthedocs.io/):

- geomeTRIC runs in Python and asks for a gradient at every step. Each one is a full
  ``gradient`` calculation of the program in ``subprograms``, with that program's own
  model and keywords, so any program with gradients works: CREST, xtb, Psi4, ORCA,
  TeraChem, ... Swap the ``ProgramSpec`` below for one of the commented alternatives.
- ``keywords`` are geomeTRIC's optimizer parameters (``geometric.params.OptParams``,
  the long ``geometric-optimize --help`` names). geomeTRIC ignores names it does not
  know, so check spelling against its documentation.
- ``coordsys``: "tric" (default, translation-rotation internal coordinates, good for
  several molecules), "dlc", "prim", or "cart".
- ``constraints`` use geomeTRIC's dictionary format with 0-based atom indices,
  distances in bohr, and angles in degrees. Scans are not supported.
- Each step in ``results.frames`` is a frame: its structure and the child's
  results (energy, gradient, charges, ...), whichever program it is.
- QM/MM: this is the way to optimize a QM region in point charges with Psi4, since
  Psi4 1.11's own optimizer cannot (see the Psi4 child below). The charges stay fixed.
"""

from qccompute import compute
from qcdata import CalcType, Model, ProgramInput, ProgramSpec, 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,
)

# The program that computes energies and gradients: CREST with GFN2-xTB.
gradient_program = ProgramSpec(
    program="crest",
    calctype=CalcType.gradient,
    model=Model(method="gfn2"),
    keywords={"threads": 1, "calculation": {"level": [{"alpb": "water"}]}},
)

# xtb instead (examples/xtb/single_point.py lists its keywords):
# gradient_program = ProgramSpec(
#     program="xtb",
#     calctype=CalcType.gradient,
#     model=Model(method="gfn2"),
#     keywords={"alpb": "water"},
# )

# Psi4 instead, here with QM/MM point charges [charge, [x, y, z]] in bohr:
# gradient_program = ProgramSpec(
#     program="psi4",
#     calctype=CalcType.gradient,
#     model=Model(method="wb97x-d", basis="def2-svp"),
#     keywords={
#         "external_potentials": [[-0.2, [0.0, 0.0, 8.0]], [0.2, [2.0, 2.0, 9.0]]],
#     },
#     cmdline_args=["-n", "1", "--memory", "2GB"],
# )

# TeraChem instead (GPU); pass propagate_wfn=True to compute() below:
# gradient_program = ProgramSpec(
#     program="terachem",
#     calctype=CalcType.gradient,
#     model=Model(method="wb97xd3", basis="def2-svp"),
#     keywords={"dftd": "d3"},
# )

prog_input = ProgramInput(
    program="geometric",
    calctype=CalcType.optimization,
    structure=structure,
    subprograms=[gradient_program],
    keywords={
        "coordsys": "tric",  # or "dlc", "prim", "cart"
        "maxiter": 100,  # maximum optimization steps
        "convergence_set": "gau",  # or "gau_loose", "gau_tight", "gau_verytight", ...
        # Individual criteria override the set (Hartree, Hartree/bohr, angstrom):
        # "convergence_energy": 1e-6,
        # "convergence_grms": 3e-4,
        # "convergence_gmax": 4.5e-4,
        # "convergence_drms": 1.2e-3,
        # "convergence_dmax": 1.8e-3,
        # "trust": 0.1,  # starting trust radius (angstrom)
        # "tmax": 0.3,  # largest trust radius (angstrom)
        #
        # Constraints (0-based atom indices; distances in bohr, angles in degrees):
        # "constraints": {
        #     "freeze": [{"type": "xyz", "indices": [0]}],  # fix the oxygen in place
        #     "set": [{"type": "angle", "indices": [1, 0, 2], "value": 104.5}],
        # },
    },
)

output = compute(
    prog_input,
    # propagate_wfn=True,  # TeraChem: start each step from the previous orbitals
    # timeout=3600,  # seconds before the whole optimization is stopped
    # collect_files=True,  # keep each child's output files in its results.files
    # executables={"crest": "/path/to/crest"},  # the child's executable
)

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

Transition state

examples/geometric/transition_state.py

"""geomeTRIC transition-state search driving another program's gradients.

Requires: ``pip install qccompute qccodec[geometric]`` (geomeTRIC itself) and the
gradient program on PATH. This example uses CREST for fast GFN2-xTB gradients:
https://github.com/crest-lab/crest/releases (official static build).

Notes for geomeTRIC transition-state searches (https://geometric.readthedocs.io/):

- ``CalcType.transition_state`` turns on geomeTRIC's saddle-point optimizer.
  geomeTRIC needs a starting Hessian, which it builds by finite differences of child
  gradients (two per Cartesian coordinate). Those displaced gradients are not
  optimizer steps, so ``results.frames`` leaves them out: it holds one frame per
  step, and each Hessian geomeTRIC computes is attached to the frame at its
  geometry (``frame.results.hessian``).
- ``"hessian"`` chooses when: ``"first"`` (the default) at the start,
  ``"first+last"`` also at the end, which is how this example checks its result,
  ``"each"`` at every step (costly).
- Start close to the saddle point. This example starts from nearly planar ammonia;
  its inversion transition state is planar, with one imaginary frequency.
- A transition state has exactly one negative Hessian eigenvalue (one imaginary
  frequency). The final frame's Hessian shows it, whatever the child program.
- The child program, coordinate systems, constraints, and QM/MM point charges work
  as in optimization.py; swap the ``ProgramSpec`` for xtb, Psi4, TeraChem, ...
- Transition-state steps use a small trust radius by default (``trust`` 0.01,
  ``tmax`` 0.03 angstrom), so expect more steps than for a minimum.
"""

import numpy as np
from qccompute import compute
from qcdata import CalcType, Model, ProgramInput, ProgramSpec, Structure

# Ammonia near its planar inversion saddle point, coordinates in bohr.
structure = Structure(
    symbols=["N", "H", "H", "H"],
    geometry=[
        [0.0, 0.0, 0.01],
        [1.88, 0.0, 0.01],
        [-0.94, 1.62, 0.01],
        [-0.94, -1.62, 0.01],
    ],
    charge=0,
    multiplicity=1,
)

# The program that computes energies and gradients: CREST with GFN2-xTB.
gradient_program = ProgramSpec(
    program="crest",
    calctype=CalcType.gradient,
    model=Model(method="gfn2"),
    keywords={"threads": 1},  # add "calculation": {"level": [{"alpb": "water"}]}
)

# Psi4 instead (any gradient program works; see optimization.py for more):
# gradient_program = ProgramSpec(
#     program="psi4",
#     calctype=CalcType.gradient,
#     model=Model(method="wb97x-d", basis="def2-svp"),
#     cmdline_args=["-n", "1", "--memory", "2GB"],
# )

prog_input = ProgramInput(
    program="geometric",
    calctype=CalcType.transition_state,
    structure=structure,
    subprograms=[gradient_program],
    keywords={
        "coordsys": "tric",  # or "dlc", "prim", "cart"
        "maxiter": 100,
        "convergence_set": "gau",  # or "gau_tight", ...
        "hessian": "first+last",  # exact Hessians at the start and the end
        # "trust": 0.01,  # starting trust radius (angstrom)
        # "tmax": 0.03,  # largest trust radius (angstrom)
    },
)

output = compute(
    prog_input,
    # propagate_wfn=True,  # TeraChem: start each step from the previous orbitals
    # timeout=3600,
)

results = output.results
print("Steps:", len(results.frames))
print("Transition-state energy (Hartree):", results.energies[-1])
print("Transition-state structure (bohr):\n", results.final_structure.geometry)

# Confirm the saddle point: exactly one negative eigenvalue of the final Hessian.
final_hessian = results.frames[-1].results.hessian
eigenvalues = np.linalg.eigvalsh(final_hessian)
print("Negative Hessian eigenvalues:", int((eigenvalues < -1e-3).sum()), "(1 for a TS)")

Many calculations at once

examples/geometric/many_calculations.py

"""geomeTRIC for many structures at once: optimizations run side by side.

Requires: what examples/geometric/optimization.py requires (geomeTRIC and CREST), and
qccompute 0.21.0 or newer.

How lists of Python calculations work (geomeTRIC, tblite):

- geomeTRIC runs in Python, and a Python calculation takes over its process's output
  streams and working directory while it runs. So when a list holds two or more,
  qccompute runs them side by side in worker processes, started with
  ``python -m qccompute._worker``: this script is never imported again and needs no
  ``if __name__ == "__main__":`` guard.
- ``max_workers`` is the budget of CPU threads, by default the cores this process may
  use. An optimization counts the threads of its gradient program (CREST ``threads``
  here). Each worker starts with ``OMP_NUM_THREADS``, ``OPENBLAS_NUM_THREADS``, and
  ``MKL_NUM_THREADS`` set to its share (unless you set them), so numpy in each worker
  does not start a thread per core.
- The gradient calculations run inside each optimization as usual.
- Failures are reported as for any list: see examples/xtb/many_calculations.py.
"""

import random
import time

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

random.seed(0)
start_geometry = [[0.0, 0.0, 0.0], [0.524, 1.687, 0.48], [1.146, -0.45, -1.354]]


def jiggled() -> Structure:
    """The starting water with every coordinate moved by up to 0.2 bohr."""
    return Structure(
        symbols=["O", "H", "H"],
        geometry=[
            [x + random.uniform(-0.2, 0.2) for x in atom] for atom in start_geometry
        ],
    )


gradient_program = ProgramSpec(
    program="crest",
    calctype=CalcType.gradient,
    model=Model(method="gfn2"),
    keywords={"threads": 1},  # threads for each gradient, counted against max_workers
)
inputs = [
    ProgramInput(
        program="geometric",
        calctype=CalcType.optimization,
        structure=jiggled(),
        subprograms=[gradient_program],
        keywords={"maxiter": 100},
    )
    for _ in range(4)
]

start = time.perf_counter()
outputs = compute(inputs)  # four workers, one optimization each
# outputs = compute(inputs, max_workers=2)  # at most two at a time
elapsed = time.perf_counter() - start

for i, output in enumerate(outputs):
    results = output.results
    print(
        f"Optimization {i}: {len(results.frames)} steps, final energy "
        f"{results.energies[-1]:.8f} Hartree"
    )
print(f"Four optimizations side by side: {elapsed:.2f} s")