CREST

CREST: write input.toml, run crest input.toml, and parse its output files.

Inputs map to CREST 3's TOML input as follows:

• keywords are the TOML document itself, so CREST's own input-file keys apply: top-level settings such as threads or runtype, and sections such as {"calculation": {"level": [{"alpb": "water"}]}} or {"cregen": {...}}. CREST stops on keys it does not know.
• model.method (gfn2, gfn1, gfn0, gfnff, ...) and the structure's charge and multiplicity are written into every [[calculation.level]]; setting charge or uhf in keywords is an error.
• runtype defaults from the calculation type (singlepoint, numhess, optimize, imtd-gc); a given runtype must match it, e.g. nci or entropy for conformer searches.
• threads defaults to the number of CPUs, at most 16. CREST itself always runs with OMP_NUM_THREADS=1 and parallelizes over its own threads.

Calculation types: energy, gradient, hessian (numerical: also frequencies and normal modes), optimization, and conformer_search. Rotamers are parsed only with the decode option {"collect_rotamers": True}; there can be thousands. CREST optimizes with the center of mass at the origin; the codec shifts optimization trajectories back into the input's frame. No input field carries atom indices.

CREST computes GFN-xTB and GFN-FF with its built-in tblite and gfnff libraries, so it needs no separate xtb. Use the official static build: conda-forge's CREST 3.0.2 crashes in its optimizer.

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

"""CREST single point: GFN-xTB or GFN-FF energy, gradient, or Hessian.

Requires: ``pip install qccompute`` (it installs qccodec) and CREST 3 on PATH, the
official static build from https://github.com/crest-lab/crest/releases (conda-forge's
CREST 3.0.2 crashes in its optimizer). CREST computes GFN-xTB and GFN-FF with its
built-in tblite and gfnff libraries, so no separate xtb is needed.

How the codec talks to CREST (https://crest-lab.github.io/crest-docs/page/documentation/inputfiles.html):

- ``keywords`` are CREST's TOML input file as a dictionary: top-level settings such
  as ``threads``, and sections such as ``calculation``. CREST stops on any key it
  does not know, so typos fail loudly.
- ``model.method`` ("gfn2", "gfn1", "gfn0", "gfnff") and the structure's charge and
  multiplicity are written into every ``[[calculation.level]]``; level settings
  (solvent, electronic temperature, accuracy) go in the one level shown below. Do
  not set ``charge`` or ``uhf`` in keywords.
- ``runtype`` comes from the calculation type ("singlepoint" or "numhess").
- The Hessian is numerical and comes with harmonic frequencies and normal modes.
- For xtb's full set of options (and g-xTB), use program "xtb" (examples/xtb/).
"""

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: CREST gets multiplicity - 1 unpaired electrons
)

prog_input = ProgramInput(
    program="crest",
    calctype=calctype,
    structure=structure,
    model=Model(method="gfn2"),  # or "gfn1", "gfn0", "gfnff"
    keywords={
        # Top-level CREST settings:
        "threads": 1,  # default: number of CPUs, at most 16
        # Settings for the method itself, one [[calculation.level]] table:
        "calculation": {
            "level": [
                {
                    "alpb": "water",  # implicit solvent; or "gbsa": "water"
                    # "etemp": 300.0,  # electronic temperature (K)
                    # "accuracy": 1.0,  # SCC accuracy; lower is tighter
                    # "maxscc": 300,  # maximum SCC iterations
                }
            ]
        },
    },
)

output = compute(
    prog_input,
    # timeout=600,  # seconds before CREST is killed
    # collect_files=True,  # keep every output file in output.results.files
    # executables={"crest": "/path/to/crest"},  # or QCCOMPUTE_CREST_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("CREST version:", results.provenance.program_version)

Optimization

examples/crest/optimization.py

"""CREST geometry optimization with its built-in ANCopt optimizer.

Requires: ``pip install qccompute`` (it installs qccodec) and CREST 3 on PATH, the
official static build from https://github.com/crest-lab/crest/releases (conda-forge's
CREST 3.0.2 crashes in its optimizer). No separate xtb is needed.

Notes for CREST optimizations (https://crest-lab.github.io/crest-docs/):

- Optimizer settings live in the ``calculation`` section of CREST's TOML input:
  ``optlev`` ("crude", "sloppy", "loose", "normal", "tight", "vtight", "extreme")
  and ``maxcycle``. Method settings (solvent, electronic temperature, ...) go in its
  ``level`` table, as in single_point.py. CREST stops on keys it does not know.
- Each step in ``results.trajectory`` has its structure and energy; the last step
  also has the gradient.
- For transition states, or an optimizer with more control, drive CREST 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="crest",
    calctype=CalcType.optimization,
    structure=structure,
    model=Model(method="gfn2"),  # or "gfn1", "gfnff"
    keywords={
        "threads": 1,  # default: number of CPUs, at most 16
        "calculation": {
            # ANCopt settings:
            "optlev": "tight",  # convergence level
            "maxcycle": 100,  # maximum optimization steps
            # The method's settings:
            "level": [
                {
                    "alpb": "water",  # implicit solvent; or "gbsa": "water"
                    # "etemp": 300.0,  # electronic temperature (K)
                }
            ],
        },
    },
)

output = compute(
    prog_input,
    # timeout=600,  # seconds before CREST is killed
    # collect_files=True,  # keep every output file (crestopt.xyz, ...) in results.files
    # executables={"crest": "/path/to/crest"},  # or QCCOMPUTE_CREST_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)
# Save the trajectory, e.g. to view in a molecule viewer:
# results.save("optimization.json"); results.final_structure.save("final.xyz")

Conformer search

examples/crest/conformer_search.py

"""CREST conformer search (iMTD-GC: metadynamics, then optimization and sorting).

Requires: ``pip install qccompute`` (it installs qccodec) and CREST 3 on PATH, the
official static build from https://github.com/crest-lab/crest/releases. No separate
xtb is needed.

Notes for CREST conformer searches (https://crest-lab.github.io/crest-docs/):

- The default workflow is iMTD-GC (``runtype = "imtd-gc"``). Other conformer
  workflows: "imtd-smtd", "entropy", "nci", "nci-mtd"; any other runtype is an error
  for this calculation type.
- Searches are long: minutes to hours for drug-sized molecules. GFN-FF ("gfnff") is
  much faster than GFN2-xTB and is a common choice for the search itself.
- ``model.method`` applies to every level of the calculation, so a GFN-FF search
  with GFN2-xTB re-ranking needs two runs: search here, then single points or
  optimizations of the conformers you keep.
- ``cregen`` settings control how structures are sorted and deduplicated: ``ewin``
  (energy window, kcal/mol), ``ethr`` (energy threshold, kcal/mol), ``rthr`` (RMSD
  threshold, angstrom), ``bthr`` (rotational-constant threshold).
- ``results.conformers`` are sorted by energy, lowest first; ``conformer_energies``
  are in Hartree. Rotamers (often thousands) are read only when asked for, with
  ``options={"collect_rotamers": True}``.
- Water has a single conformer; replace it with your molecule (Structure.open reads
  an xyz file in angstrom) and expect a longer run.
"""

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="crest",
    calctype=CalcType.conformer_search,
    structure=structure,
    model=Model(method="gfnff"),  # or "gfn2" (slower)
    keywords={
        "threads": 4,  # default: number of CPUs, at most 16
        # "runtype": "nci",  # for non-covalent clusters; default "imtd-gc"
        # "topo": False,  # allow changes in bonding topology between conformers
        "calculation": {
            "level": [
                {
                    "alpb": "water",  # implicit solvent; or "gbsa": "water"
                }
            ],
        },
        # Sorting and deduplication of the ensemble:
        "cregen": {
            "ewin": 6.0,  # keep conformers within 6 kcal/mol of the lowest
            # "ethr": 0.05,  # energy threshold for duplicates (kcal/mol)
            # "rthr": 0.125,  # RMSD threshold for duplicates (angstrom)
            # "bthr": 0.01,  # rotational-constant threshold for duplicates
        },
    },
)

output = compute(
    prog_input,
    # options={"collect_rotamers": True},  # also read crest_rotamers.xyz
    # timeout=3600,  # seconds before CREST is killed
    # collect_files=True,  # keep every output file in output.results.files
    # executables={"crest": "/path/to/crest"},  # or QCCOMPUTE_CREST_EXECUTABLE
)

results = output.results
print("Conformers:", len(results.conformers))
lowest = results.conformer_energies[0]
for i, energy in enumerate(results.conformer_energies):
    print(f"  {i}: {energy:.8f} Hartree, {(energy - lowest) * 627.5095:.2f} kcal/mol")
print("Lowest conformer (bohr):\n", results.conformers[0].geometry)
# Save the ensemble, e.g. one xyz file per conformer:
# for i, conformer in enumerate(results.conformers):
#     conformer.save(f"conformer_{i}.xyz")