"""The ORCA input writer (roadmap WP 2.2).
Translates a job specification into one ORCA input file through the keyword
model in :mod:`vibesketch.inputs.orca_keywords`, and refuses, with the
reasons, whatever that model cannot say: an unknown method or basis set, a
dispersion correction for a method that has one built in, an auxiliary basis
the approximation cannot use, frequencies for a method ORCA has no analytic
Hessian for, an optimisation with coupled cluster, a species without
electrons. It never picks a keyword the specification did not ask for.
The layout follows the qc-input-library ORCA inputs (its MANIFEST.md): the
simple-input line first, with ``TightSCF`` as the library always sets it and
``UKS`` or ``UHF`` spelt out for open shells; comment lines naming the
species and the job; ``%maxcore`` and ``%pal`` when the job names resources;
the method's blocks; the geometry in angstrom with ten decimals. No blank
lines. SMD goes through ``%cpcm``, the form ORCA 5 and 6 both read. File names
follow the library's rule, ``<compound>_<method>_<basis>[_<options>]``, with
the species id as the compound.
"""
from __future__ import annotations
import re
from collections.abc import Sequence
from typing import ClassVar
from vibesketch.document.model import Constraint, JobSpec, Species, Task
from vibesketch.inputs import orca_keywords as keywords
from vibesketch.inputs.base import InputFile, InputWriter, first_conformer, needs_conformer
__all__ = ["OrcaWriter"]
#: ORCA's simple-input keywords for each task; TDDFT and NMR also add blocks.
_TASK_KEYWORDS: dict[Task, list[str]] = {
"sp": [], # ORCA's default
"opt": ["Opt"],
"freq": ["Freq"],
"opt+freq": ["Opt", "Freq"],
"tddft": [],
"nmr": ["NMR"],
}
#: ORCA's %geom constraint letters.
_CONSTRAINT_LETTERS = {"distance": "B", "angle": "A", "dihedral": "D", "position": "C"}
#: Methods that need a /C auxiliary basis for their correlation part.
_RI_CORRELATION = {"RI-MP2", "DLPNO-CCSD", "DLPNO-CCSD(T)"}
[docs]
class OrcaWriter(InputWriter):
name: ClassVar[str] = "orca"
description: ClassVar[str] = "ORCA 6 input: DFT, HF, MP2 and coupled cluster; SP to NMR"
tasks: ClassVar[frozenset[Task]] = frozenset(_TASK_KEYWORDS)
fields: ClassVar[frozenset[str]] = frozenset(
{
"method",
"basis",
"aux_basis",
"approximation",
"dispersion",
"solvent",
"states",
"constraints",
"resources",
}
)
def _refusals(self, species: Species, spec: JobSpec) -> list[str]:
reasons = needs_conformer(self.name, species)
if species.conformers and first_conformer(species).electron_count(species.charge) == 0:
reasons.append(
f"species {species.id} has no electrons (a bare proton): there is nothing for "
"ORCA to compute; its electronic energy is zero"
)
if spec.method is None:
return reasons # the base class refuses tasks without a method
method = keywords.method_for(spec.method)
if method is None:
reasons.append(
f"orca has no method {spec.method!r} in vibe-sketch's keyword table; "
f"known: {', '.join(keywords.known_methods())}"
)
return reasons
reasons.extend(self._basis_refusals(method, spec))
reasons.extend(self._task_refusals(method, spec))
if spec.dispersion is not None:
if keywords.dispersion_keyword(spec.dispersion) is None:
reasons.append(
f"orca has no dispersion correction {spec.dispersion!r} in vibe-sketch's "
"table; known: D2, D3ZERO, D3BJ, D4"
)
elif method.has_dispersion:
reasons.append(
f"{method.keyword} includes its own dispersion correction; leave "
f"dispersion={spec.dispersion!r} out"
)
elif method.kind in ("mp2", "coupled-cluster"):
reasons.append(f"a dispersion correction is for DFT and HF, not {method.keyword}")
if spec.solvent is not None and keywords.solvent_name(spec.solvent.name) is None:
reasons.append(
f"orca's {spec.solvent.model} has no solvent {spec.solvent.name!r} in "
"vibe-sketch's table"
)
if spec.constraints and species.conformers:
atoms = len(first_conformer(species).symbols)
for constraint in spec.constraints:
beyond = [index for index in constraint.atoms if index >= atoms]
if beyond:
reasons.append(
f"a {constraint.kind} constraint names atom {beyond[0]}, but species "
f"{species.id} has atoms 0 to {atoms - 1}"
)
return reasons
def _basis_refusals(self, method: keywords.Method, spec: JobSpec) -> list[str]:
reasons: list[str] = []
if method.kind == "composite":
for field in ("basis", "aux_basis", "approximation"):
value = getattr(spec, field)
if value is not None:
reasons.append(
f"{method.keyword} brings its own basis set and approximations; "
f"leave {field}={value!r} out"
)
return reasons
if spec.basis is None:
reasons.append(f"{method.keyword} needs a basis set, for example def2-TZVP")
elif keywords.basis_family(spec.basis) is None:
reasons.append(
f"orca has no basis set {spec.basis!r} in vibe-sketch's table (def2, Pople, "
"Dunning, Jensen and x2c families)"
)
aux_kind = None
if spec.aux_basis is not None:
aux_kind = keywords.aux_basis_kind(spec.aux_basis)
if aux_kind is None:
reasons.append(
f"orca has no auxiliary basis set {spec.aux_basis!r} in vibe-sketch's table "
"(def2/J, def2/JK, def2-TZVP/C, cc-pVTZ/C, AutoAux ...)"
)
if method.keyword in _RI_CORRELATION:
if spec.approximation is not None:
reasons.append(
f"{method.keyword} with {spec.approximation} needs a /J and a /C auxiliary "
"basis set, and the specification has room for one; leave the approximation out"
)
if spec.aux_basis is None:
reasons.append(
f"{method.keyword} needs a /C auxiliary basis set, for example def2-TZVP/C"
)
elif aux_kind not in (None, "C", "auto"):
reasons.append(
f"{method.keyword} needs a /C auxiliary basis set, not {spec.aux_basis!r}"
)
return reasons
if aux_kind == "C":
reasons.append(
f"{spec.aux_basis!r} is a /C basis, for RI-MP2 and DLPNO; {method.keyword} "
"does not use one"
)
if spec.approximation is not None:
approximation = keywords.approximation_keyword(spec.approximation)
if approximation is None:
known = ", ".join(keywords.APPROXIMATIONS)
reasons.append(f"orca has no approximation {spec.approximation!r}; known: {known}")
elif method.kind in ("mp2", "coupled-cluster"):
reasons.append(
f"vibe-sketch writes {approximation} for HF and DFT, not {method.keyword}"
)
elif spec.aux_basis is None:
kinds = keywords.APPROXIMATIONS[approximation]
example = "def2/JK" if "JK" in kinds else "def2/J"
reasons.append(
f"{approximation} needs its auxiliary basis set, for example {example}"
)
elif aux_kind is not None and aux_kind not in keywords.APPROXIMATIONS[approximation]:
reasons.append(
f"{approximation} cannot use the auxiliary basis set {spec.aux_basis!r}"
)
return reasons
def _task_refusals(self, method: keywords.Method, spec: JobSpec) -> list[str]:
reasons: list[str] = []
if spec.task in ("opt", "opt+freq") and method.kind == "coupled-cluster":
reasons.append(
f"vibe-sketch does not write optimisations with {method.keyword}: optimise with "
"DFT, then run it as a single point at that geometry"
)
if spec.task in ("freq", "opt+freq") and not method.analytic_hessian:
reasons.append(
f"orca has no analytic Hessian for {method.keyword}, and vibe-sketch writes "
"Freq only where it has; compute the frequencies with HF or DFT"
)
if spec.task == "tddft" and method.kind in ("mp2", "coupled-cluster"):
reasons.append(
f"TDDFT needs an HF or DFT method, not {method.keyword}; ORCA's EOM and STEOM "
"inputs are not written yet"
)
if spec.task == "nmr" and method.kind == "coupled-cluster":
reasons.append(f"orca computes no NMR shieldings with {method.keyword}")
return reasons
[docs]
def keyword_line(self, species: Species, spec: JobSpec) -> str:
"""The simple-input line, in the library's order: reference, method,
dispersion, basis, auxiliary basis, approximation, TightSCF, task, solvent."""
method = keywords.method_for(spec.method or "")
assert method is not None # refusals ran first
reference = _reference(method, species.multiplicity)
words = [reference] if reference else []
if not (method.kind == "hf" and reference): # "UHF" alone is the method
words.append(method.keyword)
if spec.dispersion is not None:
words.append(keywords.dispersion_keyword(spec.dispersion) or spec.dispersion)
if spec.basis is not None:
words.append(spec.basis)
if spec.aux_basis is not None:
words.append(spec.aux_basis)
if spec.approximation is not None:
words.append(keywords.approximation_keyword(spec.approximation) or spec.approximation)
words.append("TightSCF")
words.extend(_TASK_KEYWORDS[spec.task])
if spec.solvent is not None:
words.append(f"CPCM({keywords.solvent_name(spec.solvent.name)})")
return "! " + " ".join(words)
[docs]
def file_name(self, species: Species, spec: JobSpec) -> str:
"""``<species>_<method>_<basis>[_<options>].inp``, by qc-input-library's naming rules."""
method = keywords.method_for(spec.method or "")
assert method is not None
reference = _reference(method, species.multiplicity)
name = method.keyword
if reference: # the library's uhf_ and uks-b3lyp_
name = reference if method.kind == "hf" else f"{reference}-{name}"
parts = [species.id, _slug(name)]
if spec.basis is not None:
parts.append(_slug(spec.basis))
tags = []
if spec.dispersion is not None:
tags.append(keywords.dispersion_keyword(spec.dispersion) or spec.dispersion)
if spec.approximation is not None:
tags.append(keywords.approximation_keyword(spec.approximation) or spec.approximation)
if spec.solvent is not None:
tags.append(f"{spec.solvent.model}-{spec.solvent.name}")
if spec.task != "sp":
tags.extend(spec.task.split("+"))
if spec.constraints:
tags.append("constrained")
parts.extend(sorted(_slug(tag) for tag in tags))
return "_".join(parts) + ".inp"
def _render(self, species: Species, spec: JobSpec) -> Sequence[InputFile]:
from vibesketch import __version__
conformer = first_conformer(species)
method = keywords.method_for(spec.method or "")
assert method is not None
what = species.name.value if species.name is not None else species.smiles or species.id
model = " / ".join(x for x in (method.keyword, spec.basis) if x)
lines = [
self.keyword_line(species, spec),
f"# {_formula(species)} / {model} — {_purpose(spec)}",
f"# {species.id}: {what}; vibe-sketch {__version__}, geometry: {conformer.method}",
*_blocks(spec),
f"* xyz {species.charge} {species.multiplicity}",
]
for symbol, (x, y, z) in zip(conformer.symbols, conformer.xyz_angstrom, strict=True):
lines.append(f"{symbol:<2}{x:14.10f}{y:15.10f}{z:15.10f}")
lines.append("*")
return [InputFile(self.file_name(species, spec), "\n".join(lines) + "\n", "text/plain")]
def _reference(method: keywords.Method, multiplicity: int) -> str | None:
"""``UKS`` or ``UHF`` for an open shell, which the library spells out; else ``None``."""
if multiplicity == 1:
return None
return "UKS" if method.kind in ("dft", "double-hybrid", "composite") else "UHF"
def _slug(text: str) -> str:
"""Lower case, every run of other characters a hyphen: def2-TZVP -> def2-tzvp."""
return re.sub(r"[^a-z0-9]+", "-", text.lower()).strip("-")
def _formula(species: Species) -> str:
if species.smiles is None:
return species.id
from rdkit.Chem import rdMolDescriptors
from vibesketch.chem.species import mol_from_smiles
return str(rdMolDescriptors.CalcMolFormula(mol_from_smiles(species.smiles)))
def _purpose(spec: JobSpec) -> str:
return {
"sp": "single point",
"opt": "geometry optimisation",
"freq": "harmonic frequencies",
"opt+freq": "optimisation, then harmonic frequencies",
"tddft": f"TDDFT, {spec.states} excited states",
"nmr": "NMR shieldings",
}[spec.task]
def _constraint_line(constraint: Constraint) -> str:
parts = [_CONSTRAINT_LETTERS[constraint.kind], *(str(i) for i in constraint.atoms)]
if constraint.value is not None:
parts.append(f"{constraint.value:.4f}")
return " { " + " ".join(parts) + " C }"
def _blocks(spec: JobSpec) -> list[str]:
"""The % blocks in the library's order: %maxcore, %pal, then the method's blocks."""
lines: list[str] = []
if spec.resources is not None:
if spec.resources.maxcore_mb is not None:
lines.append(f"%maxcore {spec.resources.maxcore_mb}")
lines.append(f"%pal nprocs {spec.resources.nprocs} end")
if spec.solvent is not None and spec.solvent.model == "SMD":
# The form ORCA 5 and 6 both read; ORCA 6 alone also takes ! SMD(solvent).
solvent = keywords.solvent_name(spec.solvent.name)
lines.extend(["%cpcm", " smd true", f' SMDsolvent "{solvent}"', "end"])
if spec.task == "tddft":
# TDA is ORCA's default; the task asks for full TDDFT.
lines.extend(["%tddft", f" nroots {spec.states}", " tda false", "end"])
if spec.constraints:
lines.extend(["%geom", " Constraints"])
lines.extend(_constraint_line(c) for c in spec.constraints)
lines.extend([" End", "end"])
return lines