vibe-sketcherDocumentation
Source code for vibesketcher.inputs.vibeqc_periodic
"""vibe-qc inputs for periodic structures: PBC, cyclic clusters and periodic MACE.
The writer turns one `PeriodicStructure` and a `PeriodicJobSpec` into a Python
script for ``vibeqc run`` (or plain ``python``). It fails closed: a request the
engine cannot run is refused with a reason, never mapped to a nearby keyword.
Models:
- ``pbc``: ``run_periodic_job`` with HF or a functional on BIPOLE, GDF or the
automatic J/K route, and a Gamma, Gamma-centred or Monkhorst-Pack mesh
(H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976)).
- ``seccm``: the semiempirical (MSINDO) cyclic cluster model; the structure's
atoms are the cluster and its periodic vectors are the CCM translations, so
1, 2 or 3 periodic vectors give CCM1D, CCM2D or CCM3D (T. Bredow,
G. Geudtner and K. Jug, J. Comput. Chem. 22, 861 (2001)).
- ``aiccm``: the ab initio cyclic cluster model, ``method="aiccm"`` with an
explicit variant and Born-von Karman torus; vibe-qc admits dim = 3 only.
- ``mace``: a periodic MACE single point or optimisation.
Any model can be preceded by a periodic MACE pre-optimisation of the atomic
positions (fixed sites honoured) or, for 3-D cells without fixed atoms, the
cell. SECCM and AICCM are experimental in vibe-qc and the input says so.
"""
from __future__ import annotations
import re
from typing import Literal
from pydantic import BaseModel, ConfigDict, Field
from vibesketcher.document.model import Document, PeriodicStructure
from vibesketcher.elements import ATOMIC_NUMBER
from vibesketcher.inputs.base import InputBundle, InputFile, WriterRefusal
from vibesketcher.units import ANGSTROM_TO_BOHR
_FUNCTIONAL = re.compile(r"^[A-Za-z0-9][A-Za-z0-9_+\-()]{0,39}$")
_BASIS = re.compile(r"^[A-Za-z0-9][A-Za-z0-9_+\-*()]{0,39}$")
_MODEL = re.compile(r"^[A-Za-z0-9][A-Za-z0-9_.\-]{0,63}$")
_AICCM_VARIANTS = ("real-gamma", "neutral-bloch", "four-center", "chi")
class _Spec(BaseModel):
model_config = ConfigDict(extra="forbid")
[docs]
class KMesh(_Spec):
"""Brillouin-zone sampling; nonperiodic directions must stay at 1."""
kind: Literal["gamma", "gamma-centred", "monkhorst-pack"] = "gamma"
mesh: tuple[int, int, int] = (1, 1, 1)
[docs]
class MaceStage(_Spec):
"""Periodic MACE relaxation run before the main job."""
relax: Literal["positions", "cell"] = "positions"
model: str = "medium-mpa-0"
fmax: float = Field(default=0.05, gt=0)
max_steps: int = Field(default=200, ge=1, le=100000)
fix_builder_sites: bool = True
fixed_sites: list[str] = Field(default_factory=list)
device: Literal["cpu", "cuda", "mps"] = "cpu"
accept_academic_license: bool = False
[docs]
class PeriodicJobSpec(_Spec):
"""What to compute for one periodic structure."""
model: Literal["pbc", "seccm", "aiccm", "mace"] = "pbc"
task: Literal["sp", "opt"] = "sp"
method: str = "hf"
basis: str | None = "sto-3g"
#: J/K route; None picks BIPOLE for 3-D cells and vibe-qc's automatic
#: (rigorous 2-D/1-D Coulomb) route otherwise.
jk: Literal["bipole", "gdf", "auto"] | None = None
kpoints: KMesh = Field(default_factory=KMesh)
charge: int | None = None
multiplicity: int | None = Field(default=None, ge=1)
madelung: bool = True
max_iter: int | None = Field(default=None, ge=1, le=10000)
variant: str | None = None
lattice_extension: tuple[int, int, int] | None = None
preoptimize: MaceStage | None = None
mace: MaceStage | None = None
def _refusals(structure: PeriodicStructure, spec: PeriodicJobSpec) -> list[str]:
reasons = []
dim = sum(structure.pbc)
if structure.pbc != tuple(i < dim for i in range(3)):
reasons.append("vibe-qc needs the periodic directions first (pbc T, T, F or T, F, F)")
if any(site.occupancy != 1 for site in structure.sites):
reasons.append("partially occupied sites have no vibe-qc input")
multiplicity = spec.multiplicity or structure.multiplicity or 1
electrons = sum(ATOMIC_NUMBER[s.symbol] for s in structure.sites) - (
spec.charge if spec.charge is not None else structure.charge or 0
)
if (electrons + multiplicity - 1) % 2:
reasons.append(f"{electrons} electrons cannot have multiplicity {multiplicity}")
method = spec.method.lower()
if spec.model in ("pbc", "aiccm"):
if not spec.basis or not _BASIS.match(spec.basis):
reasons.append(f"{spec.model} needs a basis set name")
if method != "hf" and not _FUNCTIONAL.match(spec.method):
reasons.append(f"cannot write functional {spec.method!r}")
if spec.model == "pbc":
if spec.jk == "bipole" and dim != 3:
reasons.append(
"BIPOLE is a bulk (dim = 3) Coulomb builder; use jk 'auto' for a true "
f"{dim}-D model or build the slab with periodicity '3d' (vacuum box)"
)
for axis in range(3):
if axis >= dim and spec.kpoints.mesh[axis] != 1:
reasons.append("k-points along a nonperiodic direction must be 1")
if spec.kpoints.kind == "gamma" and spec.kpoints.mesh != (1, 1, 1):
reasons.append("a Gamma-only calculation takes no mesh; use gamma-centred")
if any(n < 1 or n > 64 for n in spec.kpoints.mesh):
reasons.append("k-mesh entries must be between 1 and 64")
elif spec.kpoints != KMesh():
reasons.append(f"{spec.model} samples no k-points; the cluster size sets the sampling")
if spec.model == "seccm":
if method not in ("msindo", "hf"):
reasons.append("SECCM is written for the MSINDO Hamiltonian only")
if multiplicity != 1:
reasons.append("SECCM (MSINDO) is closed-shell only")
if spec.task != "sp":
reasons.append("SECCM (MSINDO) runs single points only")
if spec.model == "aiccm":
if dim != 3:
reasons.append(
"vibe-qc admits AICCM for three-dimensional cells only (vibe-qc#1089); "
"use SECCM or a 3-D vacuum box"
)
if spec.variant not in _AICCM_VARIANTS:
reasons.append(f"AICCM needs variant, one of {', '.join(_AICCM_VARIANTS)}")
if spec.lattice_extension is None or any(n < 1 for n in spec.lattice_extension):
reasons.append("AICCM needs a lattice_extension torus such as (2, 2, 2)")
if spec.task != "sp":
reasons.append("AICCM inputs are single points")
elif spec.variant is not None or spec.lattice_extension is not None:
reasons.append("variant and lattice_extension belong to AICCM")
if spec.model == "mace":
if spec.preoptimize is not None:
reasons.append("a MACE job needs no MACE pre-optimisation")
if spec.basis not in (None, "sto-3g") or method not in ("hf", "mace"):
reasons.append("a MACE job takes no method or basis")
elif spec.mace is not None:
reasons.append("mace options belong to model 'mace'; use preoptimize for a pre-stage")
for stage in (spec.preoptimize, spec.mace):
if stage is None:
continue
if not _MODEL.match(stage.model):
reasons.append(f"cannot write MACE model {stage.model!r}")
if stage.relax == "cell" and (dim != 3 or _fixed(structure, stage)):
reasons.append("MACE cell relaxation needs a 3-D cell without fixed atoms")
unknown = set(stage.fixed_sites) - {s.id for s in structure.sites}
if unknown:
reasons.append(f"fixed sites not in the structure: {', '.join(sorted(unknown))}")
return reasons
def _fixed(structure: PeriodicStructure, stage: MaceStage) -> list[int]:
ids = [s.id for s in structure.sites]
wanted = set(stage.fixed_sites)
builder = structure.provenance.get("builder")
if stage.fix_builder_sites and isinstance(builder, dict):
recorded = builder.get("fixed_sites")
if isinstance(recorded, list):
wanted.update(str(x) for x in recorded)
return [i for i, sid in enumerate(ids) if sid in wanted]
[docs]
def refusals(structure: PeriodicStructure, spec: PeriodicJobSpec) -> list[str]:
return _refusals(structure, spec)
def _stem(structure: PeriodicStructure, spec: PeriodicJobSpec) -> str:
parts = [structure.id, spec.model]
if spec.model in ("pbc", "aiccm"):
parts += [spec.method.lower(), (spec.basis or "").lower()]
if spec.model == "pbc":
parts.append("k" + "x".join(str(n) for n in spec.kpoints.mesh))
return re.sub(r"[^A-Za-z0-9_.-]+", "-", "-".join(p for p in parts if p))
def _float(x: float) -> str:
return f"{x: .10f}"
def _mace_lines(stage: MaceStage) -> list[str]:
return [
f"mace = MLIPOptions(model={stage.model!r}, device={stage.device!r}, "
f"accept_academic_license={stage.accept_academic_license!r})"
]
[docs]
def write_periodic_input(structure: PeriodicStructure, spec: PeriodicJobSpec) -> InputBundle:
"""Render the vibe-qc script for one periodic structure, or refuse."""
from vibesketcher import __version__
reasons = _refusals(structure, spec)
if reasons:
raise WriterRefusal("vibeqc-periodic", reasons)
dim = sum(structure.pbc)
stem = _stem(structure, spec)
charge = spec.charge if spec.charge is not None else structure.charge or 0
multiplicity = spec.multiplicity or structure.multiplicity or 1
method = spec.method.lower()
builder = structure.provenance.get("builder")
made = (
f"{builder.get('kind', 'structure')}" if isinstance(builder, dict) else "periodic structure"
)
status = {
"pbc": "Periodic SCF",
"seccm": "Experimental: semiempirical cyclic cluster model (MSINDO)",
"aiccm": "Experimental: ab initio cyclic cluster model",
"mace": "Periodic MACE machine-learned potential; no wavefunction or accuracy claim",
}[spec.model]
rows = [
f" ({ATOMIC_NUMBER[s.symbol]:>3}, "
+ ", ".join(_float(x) for x in structure.cartesian(s))
+ f"), # {s.id}"
for s in structure.sites
]
lattice = [
" [" + ", ".join(_float(x) for x in row) + "]," for row in structure.lattice_angstrom
]
lines = [
f'"""{structure.id}: {len(structure.sites)} sites, {dim}D {made}, '
f"{spec.model} {spec.task}.",
"",
f"Written by vibe-sketcher {__version__}. {status}.",
"Converge k-mesh, cluster size, slab thickness and vacuum before drawing conclusions.",
f"Run: vibeqc run {stem}.py",
'"""',
"",
"import vibeqc as vq",
"",
f"ANGSTROM_TO_BOHR = {ANGSTROM_TO_BOHR!r} # CODATA 2018, as vibe-qc uses",
"",
"# Lattice row vectors a1, a2, a3 in angstrom; the first "
f"{dim} {'is' if dim == 1 else 'are'} periodic.",
"LATTICE = [",
*lattice,
"]",
"# Atomic number and Cartesian position in angstrom.",
"ATOMS = [",
*rows,
"]",
"",
"",
"def bohr(xyz):",
" return [x * ANGSTROM_TO_BOHR for x in xyz]",
"",
"",
"system = vq.PeriodicSystem(",
f" {dim},",
" lattice_vectors=[bohr(row) for row in LATTICE],",
" unit_cell=[vq.Atom(z, bohr(xyz)) for z, *xyz in ATOMS],",
f" charge={charge},",
f" multiplicity={multiplicity},",
")",
"",
]
stage = spec.preoptimize
if stage is not None:
fixed = _fixed(structure, stage)
lines += [
"from vibeqc.mlip.mace import (",
" optimize_periodic_mace_cell,",
" optimize_periodic_mace_positions,",
")",
"from vibeqc.mlip.options import MLIPOptions",
"",
f"# Pre-optimise the {stage.relax} with periodic MACE before the main job.",
*_mace_lines(stage),
]
if stage.relax == "cell":
lines += [
f"relaxed = optimize_periodic_mace_cell(system, mace, fmax={stage.fmax!r}, "
f"max_steps={stage.max_steps!r})",
]
else:
lines += [
f"FIXED = {fixed!r} # site indices held fixed",
f"relaxed = optimize_periodic_mace_positions(system, mace, fmax={stage.fmax!r}, "
f"max_steps={stage.max_steps!r}, fixed_indices=FIXED)",
]
lines += ["system = relaxed.system", ""]
if spec.model == "pbc":
if spec.kpoints.kind == "gamma":
kpoints = "vq.KPoints.gamma(system)"
elif spec.kpoints.kind == "gamma-centred":
kpoints = f"vq.KPoints.gamma_centred(system, {list(spec.kpoints.mesh)!r})"
else:
kpoints = f"vq.KPoints.monkhorst_pack(system, {list(spec.kpoints.mesh)!r})"
reference = "U" if multiplicity > 1 else "R"
scf = f"{reference}HF" if method == "hf" else f"{reference}KS"
arguments = [
("method", repr(scf)),
*([("functional", repr(spec.method))] if method != "hf" else []),
("jk_method", repr(spec.jk or ("bipole" if dim == 3 else "auto"))),
("kpoints", "kpoints"),
]
if spec.task == "opt":
fixed_opt = _fixed(structure, MaceStage())
arguments.append(("optimize", "True"))
if fixed_opt:
arguments.append(("optimize_frozen_indices", repr(fixed_opt)))
if spec.max_iter:
arguments.append(("max_iter", str(spec.max_iter)))
lines += [
f"basis = vq.BasisSet(system.unit_cell_molecule(), {spec.basis!r})",
f"kpoints = {kpoints}",
"result = vq.run_periodic_job(",
" system,",
" basis,",
*(f" {k}={v}," for k, v in arguments),
f" output={stem!r},",
")",
'print(f"E = {result.energy:.10f} Ha converged = {result.converged}")',
]
elif spec.model == "aiccm":
lines += [
f"basis = vq.BasisSet(system.unit_cell_molecule(), {spec.basis!r})",
"result = vq.run_periodic_job(",
" system,",
" basis,",
' method="aiccm",',
f" variant={spec.variant!r},",
f" functional={None if method == 'hf' else spec.method!r},",
f" aiccm_lattice_extension={tuple(spec.lattice_extension or ())!r},",
*([f" max_iter={spec.max_iter},"] if spec.max_iter else []),
f" output={stem!r},",
")",
'print(f"E = {result.energy:.10f} Ha converged = {result.converged}")',
]
elif spec.model == "seccm":
translations = [list(row) for row in structure.lattice_angstrom[:dim]]
lines += [
"from vibeqc.semiempirical.methods.msindo_ccm import CCMOptions",
"",
"# The cyclic cluster is every atom of the cell; its periodic vectors are the",
f"# CCM translations in angstrom ({dim} vectors: CCM{dim}D).",
"TRANSLATIONS = [",
*(" [" + ", ".join(_float(x) for x in row) + "]," for row in translations),
"]",
"cluster = vq.Molecule(",
" [vq.Atom(z, bohr(xyz)) for z, *xyz in ATOMS],",
f" charge={charge},",
f" multiplicity={multiplicity},",
")",
"result = vq.run_job(",
" cluster,",
' method="seccm",',
" ccm_options=CCMOptions(",
" translations=TRANSLATIONS,",
f" madelung={spec.madelung!r},",
*([f" max_iter={spec.max_iter},"] if spec.max_iter else []),
" ),",
f" output={stem!r},",
")",
'print(f"E = {result.energy:.10f} Ha")',
]
else:
stage = spec.mace or MaceStage()
fixed = _fixed(structure, stage)
lines += [
"from vibeqc.mlip.options import MLIPOptions",
"",
*_mace_lines(stage),
"result = vq.run_periodic_job(",
" system,",
" None,",
' method="mace",',
" mlip_options=mace,",
]
if spec.task == "opt":
lines.append(" optimize=True,")
if stage.relax == "cell":
lines.append(" optimize_cell=True,")
elif fixed:
lines.append(f" optimize_frozen_indices={fixed!r},")
lines += [f" output={stem!r},", ")", 'print(f"E = {result.energy:.10f} Ha")']
content = "\n".join(lines) + "\n"
return InputBundle(
program="vibeqc",
species=structure.id,
files=(InputFile(f"{stem}.py", content, "text/x-python"),),
)
[docs]
def write_periodic_inputs(document: Document, structure: str, spec: PeriodicJobSpec) -> InputBundle:
"""The script for the named periodic structure of a document."""
found = {s.id: s for s in document.periodic_structures}
if structure not in found:
raise WriterRefusal("vibeqc-periodic", [f"periodic structure {structure!r} not found"])
return write_periodic_input(found[structure], spec)