TeraChem

TeraChem: write tc.in, run terachem tc.in, and parse stdout and scratch files.

Inputs map to TeraChem as follows:

• model.method and model.basis become method and basis; the structure's charge and multiplicity become charge and spinmult. Setting any of these, or run, in keywords is an error.
• Every other keyword becomes a tc.in line, key value, so TeraChem's own keyword names apply (scf, convthre, dftd, pcm, epsilon, ...). Values are written in lowercase, which suits TeraChem's booleans and options but means case-sensitive values such as file names must be lowercase. TeraChem ignores keywords it does not know, printing only "Unrecognized or unnecessary keywords" in its output, so check the logs when a setting seems to have no effect.
• cmdline_args are appended to the terachem command.

Calculation types: energy, gradient, hessian (run frequencies; the Hessian and gradient are parsed, frequencies and normal modes are not yet), and optimization, which requires new_minimizer yes in keywords. TeraChem refuses Hessians away from a stationary point unless mincheck false is set. TeraChem keeps the input's coordinates and atom order. No input field carries atom indices.

The orbital files c0 (or ca0/cb0) can seed the next calculation: qccompute's propagate_wfn=True does this between geomeTRIC steps. A run is a failure unless stdout ends with Job finished:.

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/terachem/single_point.py

"""TeraChem single point: energy, gradient, or Hessian.

Requires: ``pip install qccompute`` (it installs qccodec) and ``terachem`` on PATH.
TeraChem runs only on NVIDIA CUDA GPUs and needs a license file (from PetaChem). A
Docker launcher script that provides ``terachem`` without an install, with the license
given by TERACHEM_LICENSE_PATH, comes with the images at
https://hub.docker.com/r/mtzgroup/terachem.

How the codec talks to TeraChem (manual: https://mtzweb.stanford.edu/terachem/,
or the PDF manual shipped in the TeraChem ``doc`` directory):

- The codec writes ``tc.in``. ``run``, ``coordinates``, ``charge``, ``spinmult``,
  ``method``, and ``basis`` come from ``calctype``, ``structure``, and ``model``;
  setting any of them in ``keywords`` is an error.
- Every other ``keywords`` entry becomes one ``key value`` line in ``tc.in``, so use
  TeraChem's keyword names. Values are written lowercased (``True`` becomes ``true``);
  TeraChem itself expects ``yes``/``no`` for switches. TeraChem ignores unknown or
  misspelled keywords with only an "Unrecognized or unnecessary keywords" warning
  in ``output.logs``, so check there when a setting seems to have no effect.
- Unrestricted methods are spelled in the method name: "uhf", "ub3lyp", ...;
  restricted open shell is "rohf", "rob3lyp".
- Hessians (``run frequencies``) come from finite differences of analytic
  gradients, so they take 6N gradients; the gradient at the input geometry is
  returned as well. TeraChem stops with "Maximum component of gradient is too large"
  unless the structure is a stationary point; ``mincheck false`` turns that check
  off (set below because this example water is not at a minimum; energy and
  gradient runs list it as an unnecessary keyword, which is harmless).
- ``purify no`` is set below because density purification is TeraChem's default for
  large systems and can fail to converge for small molecules.
"""

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

prog_input = ProgramInput(
    program="terachem",
    calctype=calctype,
    structure=structure,
    model=Model(method="wb97xd3", basis="def2-svp"),  # wB97X-D3; e.g. "b3lyp", "ub3lyp"
    keywords={
        # Any tc.in keyword, exactly as in TeraChem's manual:
        "purify": "no",
        "convthre": 1e-6,  # SCF convergence (wavefunction error)
        "threall": 1e-12,  # integral screening threshold
        "maxit": 100,  # SCF iterations
        "mincheck": "false",  # allow Hessians away from a stationary point
        # "precision": "double",  # "single", "mixed" (default), "dynamic", "double"
        # "dftd": "d3",  # empirical dispersion for DFT methods
        # "gpus": "1 0",  # number of GPUs, then the device indices to use
        #
        # Implicit solvent (COSMO); "epsilon" is the dielectric constant:
        "pcm": "cosmo",
        "epsilon": 78.39,  # water
        # "pcm_scale": 1.2,  # scale factor for the cavity atomic radii
        #
        # Excited states (TDDFT/CIS); the gradient is for state "cistarget":
        # "cis": "yes",
        # "cisnumstates": 3,
        # "cistarget": 1,
    },
)

output = compute(
    prog_input,
    # timeout=3600,  # seconds before TeraChem is killed
    # collect_files=True,  # keep every output file in output.results.files
    # collect_wfn=True,  # keep the orbitals (scr*/c0) for a later guess
    # executables={"terachem": "/path/to/terachem"},  # or QCCOMPUTE_TERACHEM_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("Hessian (Hartree/bohr^2), shape:", results.hessian.shape)
print("Number of orbitals:", results.calcinfo_nmo)
if results.extras.get("excited_states"):
    print("Excited states:", results.extras["excited_states"])
print("TeraChem version:", results.provenance.program_version)

Optimization

examples/terachem/optimization.py

"""TeraChem geometry optimization with its built-in optimizer.

Requires: ``pip install qccompute`` (it installs qccodec) and ``terachem`` on PATH.
TeraChem runs only on NVIDIA CUDA GPUs and needs a license file (from PetaChem).

Notes for TeraChem optimizations (manual: https://mtzweb.stanford.edu/terachem/):

- Only TeraChem's newer optimizer is supported: ``"new_minimizer": "yes"`` is
  required, and the codec raises an error without it.
- Optimizer settings are ordinary tc.in keywords: ``nstep``, ``min_coordinates``,
  ``min_converge_e``, ``min_converge_grms``, ``min_converge_gmax``, ...
- TeraChem does not stop on a misspelled keyword: it prints "Unrecognized or
  unnecessary keywords: These are being ignored." in its output and carries on.
  Check ``output.logs`` for that warning when a setting seems to have no effect.
- Every step in ``results.trajectory`` carries its structure, energy, and gradient,
  read from TeraChem's ``optim.xyz``/``optim_trajectory.xyz`` and its output.
- The first step is your input structure: TeraChem does not move or reorient it.
- Transition-state searches are not supported by this codec; use geomeTRIC with
  TeraChem gradients (``program="geometric"``) instead.
"""

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="terachem",
    calctype=CalcType.optimization,
    structure=structure,
    model=Model(method="wb97xd3", basis="def2-svp"),
    keywords={
        "purify": "no",
        "new_minimizer": "yes",  # required
        "nstep": 50,  # maximum optimization steps
        # "min_coordinates": "cartesian",  # default: TRIC internal coordinates
        # Convergence criteria (TeraChem's defaults shown):
        # "min_converge_grms": 3.0e-4,  # RMS gradient (Hartree/bohr)
        # "min_converge_gmax": 4.5e-4,  # largest gradient component
        # "min_converge_e": 1.0e-6,  # energy change (Hartree)
        #
        # Implicit solvent (COSMO):
        # "pcm": "cosmo",
        # "epsilon": 78.39,
    },
)

output = compute(prog_input)  # timeout=..., collect_files=..., executables=...

results = output.results
print("Steps:", len(results.trajectory))
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")