Sella¶

Sella: minimizations, transition-state searches, and IRCs driving another program.

Sella (https://github.com/zadorlab/sella) is a Python optimizer built on ASE. It asks an ASE calculator for energies and forces; here that calculator is qccodec.ase.ProgramCalculator, so 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, and no ASE calculator is involved.

Inputs map to Sella as follows:

subprograms holds one ProgramSpec with calctype="gradient": the child program, its model, and its own keywords (passed to that program unchanged). Optionally a second one with calctype="hessian" becomes Sella's hessian_function: Sella then starts from that program's Hessian instead of estimating curvature from gradients, and recomputes it every diag_every_n steps if that keyword is set. The Hessian program needs analytic or cheap Hessians to pay off (xtb, AIMNet2, MACE, ...).
keywords are Sella's own, passed unchanged in Sella's units (eV and angstrom): fmax and steps go to Sella.run, everything else (eta, gamma, internal, delta0, diag_every_n, nsteps_per_diag, ...) to the Sella constructor. Sella's defaults apply to anything left out, including fmax=0.05 eV/angstrom and an effectively unlimited steps, so set steps.
constraints maps Sella's Constraints methods to their calls: fix_translation, fix_bond, fix_angle, fix_dihedral. Each call is the first argument (atom indices) or a dictionary of keyword arguments, e.g. {"fix_translation": [0], "fix_bond": [[0, 1], {"indices": [1, 2], "target": 1.0}]}. Atom indices are 0-based and targets in angstrom and degrees, as in Sella.
The codec sets order from the calculation type and keeps Sella's log on stdout; order, logfile, trajectory, restart, append_trajectory, hessian_function, and cmdline_args are not accepted.

Calculation types: optimization (order=0), transition_state (order=1), and irc (Sella's IRC from the transition state in structure, both directions). For an IRC, fmax, fmax_inner, and steps go to IRC.run and the rest (dx, keep_going, eta, gamma, ...) to IRC; direction is set by the codec. The result is a PathData: frames from the end of the reverse branch, through the transition state (ts_index), to the end of the forward branch. Two corrections to Sella keep the path meaningful. Sella's IRC keeps overall translation and rotation, so a branch drifts and can spin at its minimum: each frame is moved onto the previous one by the mass-weighted best fit, its gradient rotated with it, which keeps the path in the input's frame. And a branch ends when the curvature of the vibrations is positive, not of the full Hessian, whose near-zero rigid-motion eigenvalues noise can push negative (Sella would otherwise keep rotating the molecule until steps runs out). With constraints, Sella's own test applies.

The trajectory holds one frame per optimizer step. Gradients Sella evaluates at displaced geometries to probe the curvature are not steps and are left out, unless one fails. Each Hessian from the Hessian program is attached to the frame at its geometry (frame.results.hessian). Sella removes overall translation and rotation from its steps, so every frame stays in the input's frame. A run that uses up steps without meeting fmax fails, with its trajectory attached.

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/sella/optimization.py

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

Requires: ``pip install qccompute qccodec[sella]`` (Sella and ASE) and the gradient
program on PATH. This example uses xtb for fast GFN2-xTB gradients:
``conda install -c conda-forge xtb``.

How it works (https://github.com/zadorlab/sella):

- Sella 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: xtb, CREST, Psi4, ORCA,
  TeraChem, AIMNet2, ... ASE is only the interface Sella talks to; no ASE calculator
  runs.
- ``keywords`` are Sella's own, in its units (eV, angstrom). ``fmax`` and ``steps``
  go to ``Sella.run``; everything else to the ``Sella`` constructor. Sella's defaults
  apply to anything left out: ``fmax`` 0.05 eV/angstrom and no practical step limit,
  so set ``steps``. A run that uses up ``steps`` fails, with its trajectory attached.
- ``constraints`` are calls to Sella's ``Constraints`` methods with 0-based atom
  indices; targets in angstrom and degrees.
- For minimizations Sella estimates curvature from gradients. A Hessian program in
  ``subprograms`` is used only with ``"eig": True`` (see transition_state.py).
- Each step in ``results.frames`` is a frame: its structure and the child's results
  (energy, gradient, ...). The child's program, model, and keywords are in
  ``output.input_data.subprograms``.
"""

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: xtb with GFN2-xTB.
gradient_program = ProgramSpec(
    program="xtb",
    calctype=CalcType.gradient,
    model=Model(method="gfn2"),
    keywords={"alpb": "water"},  # xtb's own keywords; see examples/xtb/
)

prog_input = ProgramInput(
    program="sella",
    calctype=CalcType.optimization,
    structure=structure,
    subprograms=[gradient_program],
    keywords={
        "fmax": 0.01,  # largest force at convergence (eV/angstrom)
        "steps": 100,  # maximum optimizer steps
        # "internal": True,  # step in internal coordinates
        #
        # Constraints (0-based atom indices; targets in angstrom and degrees):
        # "constraints": {
        #     "fix_translation": [0],  # keep the oxygen in place
        #     "fix_angle": [{"indices": [1, 0, 2], "target": 104.5}],
        # },
    },
)

output = compute(
    prog_input,
    # propagate_wfn=True,  # TeraChem: start each step from the previous orbitals
    # executables={"xtb": "/path/to/xtb"},  # 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("Sella version:", results.provenance.program_version)

Transition state¶

examples/sella/transition_state.py

"""Sella transition-state search driving another program's gradients and Hessians.

Requires: ``pip install qccompute qccodec[sella]`` (Sella and ASE) and the child
program on PATH. This example uses xtb (GFN2-xTB): ``conda install -c conda-forge xtb``.

Notes for Sella transition-state searches (https://github.com/zadorlab/sella):

- ``CalcType.transition_state`` runs Sella's saddle-point search (``order=1``).
- Without a Hessian program, Sella estimates the curvature from extra gradients at
  displaced geometries; those are not steps and stay out of ``results.frames``.
- With a second ``ProgramSpec`` of ``calctype="hessian"`` in ``subprograms`` (here
  xtb's Hessian), Sella starts from that program's Hessian instead, which usually
  takes fewer steps. ``"diag_every_n": N`` recomputes it every N steps. Each Hessian
  is attached to the frame at its geometry (``frame.results.hessian``).
- Start close to the saddle point. This example starts from a bent HCN, between HCN
  and HNC; the isomerization transition state has one imaginary frequency.
- ``fmax`` and ``steps`` are Sella's (eV/angstrom); the defaults are Sella's.
"""

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

# H between C and N, coordinates in bohr.
structure = Structure(
    symbols=["C", "N", "H"],
    geometry=[[0.0, 0.0, 0.0], [2.2, 0.0, 0.0], [1.1, 1.9, 0.0]],
    charge=0,
    multiplicity=1,
)

model = Model(method="gfn2")
gradient_program = ProgramSpec(program="xtb", calctype=CalcType.gradient, model=model)
hessian_program = ProgramSpec(program="xtb", calctype=CalcType.hessian, model=model)

prog_input = ProgramInput(
    program="sella",
    calctype=CalcType.transition_state,
    structure=structure,
    subprograms=[gradient_program, hessian_program],  # drop the Hessian to estimate it
    keywords={
        "fmax": 0.001,  # largest force at convergence (eV/angstrom)
        "steps": 200,  # maximum optimizer steps
        # "diag_every_n": 5,  # recompute the Hessian every 5 steps
        # "internal": True,  # step in internal coordinates
    },
)

output = compute(prog_input)
results = output.results
print("Steps:", len(results.frames))
print("Energies (Hartree):", results.energies)
print("Transition state (bohr):\n", results.final_structure.geometry)

# Check: a transition state has exactly one imaginary frequency.
check = compute(
    ProgramInput(
        program="xtb",
        calctype=CalcType.hessian,
        structure=results.final_structure,
        model=model,
    )
)
print("Lowest frequencies (cm^-1):", check.results.freqs_wavenumber[:3])

IRC¶

examples/sella/irc.py

"""Sella intrinsic reaction coordinate (IRC) from a transition state.

Requires: ``pip install qccompute qccodec[sella]`` (Sella and ASE) and the child
program on PATH. This example uses xtb (GFN2-xTB): ``conda install -c conda-forge xtb``.

Notes for Sella IRCs (https://github.com/zadorlab/sella):

- ``CalcType.irc`` follows the steepest-descent path in mass-weighted coordinates
  down both sides of the transition state in ``structure``: Sella's ``IRC``, run in
  the reverse and then the forward direction from one analysis of the saddle.
- ``results`` is a ``PathData``: ``frames`` run from the end of the reverse branch,
  through the transition state (``ts_index``), to the end of the forward branch. Each
  frame holds its structure and the child's energy and gradient;
  ``forward_barrier`` and ``reverse_barrier`` are the transition state's energy above
  the first and last frames, in Hartree.
- A ``hessian`` subprogram (here xtb's Hessian) gives Sella the transition state's
  curvature instead of estimating it from gradients; it is attached to the
  transition-state frame.
- ``fmax``, ``fmax_inner``, ``steps`` go to ``IRC.run``; ``dx`` (step size,
  angstrom * amu^1/2), ``keep_going`` (continue past an inner-loop failure),
  ``eta``, ``gamma`` to ``IRC``. Sella's defaults apply to anything left out; set
  ``steps``.
- The codec keeps the path in the input's frame (Sella's IRC lets the molecule drift
  and spin) and ends a branch when the curvature of the vibrations, not the rigid
  rotations, is positive.
- An IRC stops short of the minima. Optimize the end structures before comparing
  energies (``path.frames[0].structure``, ``path.frames[-1].structure``).
"""

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

# The HCN <-> HNC transition state (GFN2-xTB), coordinates in bohr.
transition_state = Structure(
    symbols=["C", "N", "H"],
    geometry=[
        [0.017, 0.036, 0.0],
        [2.283, -0.137, 0.0],
        [1.0, 2.0, 0.0],
    ],
    charge=0,
    multiplicity=1,
)

model = Model(method="gfn2")
prog_input = ProgramInput(
    program="sella",
    calctype=CalcType.irc,
    structure=transition_state,
    subprograms=[
        ProgramSpec(program="xtb", calctype=CalcType.gradient, model=model),
        ProgramSpec(program="xtb", calctype=CalcType.hessian, model=model),
    ],
    keywords={
        "fmax": 0.01,  # force at which a branch counts as at its minimum (eV/angstrom)
        "steps": 300,  # maximum steps per direction
        # "dx": 0.1,  # step size (angstrom * amu^1/2)
        # "keep_going": True,  # step on past an inner-loop convergence failure
    },
)

output = compute(prog_input)
path = output.results
print("Frames:", len(path.frames), "transition state at", path.ts_index)
print(f"Forward barrier: {path.forward_barrier * 627.5095:.1f} kcal/mol")
print(f"Reverse barrier: {path.reverse_barrier * 627.5095:.1f} kcal/mol")
print("Reverse end (bohr):\n", path.frames[0].structure.geometry)
print("Forward end (bohr):\n", path.frames[-1].structure.geometry)
# Save the path, e.g. to view in a molecule viewer: path.save("irc.xyz")