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).

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.
- Every step in ``results.trajectory`` is the child's full output (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.trajectory))
print("Energies (Hartree):", results.energies)
print("Final gradient (Hartree/bohr):\n", results.trajectory[-1].results.gradient)
print("Final structure (bohr):\n", results.final_structure.geometry)
print("Gradient program:", results.trajectory[-1].input_data.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.trajectory`` 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.trajectory))
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.trajectory[-1].results.hessian
eigenvalues = np.linalg.eigvalsh(final_hessian)
print("Negative Hessian eigenvalues:", int((eigenvalues < -1e-3).sum()), "(1 for a TS)")