vibe-sketcherDocumentation
Source code for vibesketcher.document.solid_energetics
"""Adsorption, surface and neutral defect energies from linked periodic runs (#195).
A *plan* fixes the members of an energy difference, writes one vibe-qc input
per member with the same Hamiltonian, and records the formula, the
normalisation and every assumption. *Evaluation* reads the members' result
QVFs, checks that each result is the planned structure computed with the
planned method, and returns the energy with units. Nothing is inferred:
- adsorption: E_ads = [E(slab + n X) - E(slab) - n E(X, gas)] / n, with the
complex built from that slab and its composition differing by exactly n X;
- surface: gamma = [E(slab) - N E(bulk per formula unit)] / (2 A) for a
stoichiometric slab with identical faces (Boettger, Phys. Rev. B 49, 16798
(1994) discusses the bulk-reference choice); A is the in-plane area;
- defect formation (neutral only): E_f = E(defective) - E(pristine) -
sum_i dn_i mu_i, with every reservoir mu_i stated (a value and its source,
or a gas species and its atom count, e.g. mu_O = E(O2)/2). Charged defects
need a Fermi level, potential alignment and a finite-size correction
(Freysoldt, Neugebauer and Van de Walle, Phys. Rev. Lett. 102, 016402
(2009)); they are refused here rather than reported uncorrected.
A gas-phase molecule is computed as a molecule with the same method and basis
(``seccm`` members pair with molecular MSINDO). k-meshes, slab thickness,
vacuum and cluster size are recorded, not converged, by a plan.
"""
from __future__ import annotations
import io
import json
import math
import re
import zipfile
from collections import Counter
from typing import Any, Literal
import numpy as np
from pydantic import BaseModel, ConfigDict, Field
from vibesketcher.document import crystal
from vibesketcher.document.model import Document, PeriodicStructure, Species
from vibesketcher.elements import ATOMIC_NUMBER
from vibesketcher.inputs.base import InputBundle, InputFile, WriterRefusal
from vibesketcher.inputs.vibeqc_periodic import PeriodicJobSpec, write_periodic_input
from vibesketcher.units import ANGSTROM_TO_BOHR
PLAN_SCHEMA = "vibe-solid-energy/1"
HARTREE_EV = 27.211386245988 # CODATA 2018
EV_KJMOL = 96.48533212331
EV_PER_A2_TO_J_PER_M2 = 16.021766340
[docs]
class Reservoir(BaseModel):
"""A chemical potential: a stated value, or half (1/atoms) of a gas species' energy."""
model_config = ConfigDict(extra="forbid")
value_ev: float | None = None
source: str | None = None
species: str | None = None
atoms: int = Field(default=1, ge=1)
[docs]
class Member(BaseModel):
model_config = ConfigDict(extra="forbid")
role: str
structure: str | None = None
species: str | None = None
formula: dict[str, int]
input: str
method: str
[docs]
class SolidEnergyPlan(BaseModel):
model_config = ConfigDict(extra="forbid")
schema_id: Literal["vibe-solid-energy/1"] = Field(default="vibe-solid-energy/1", alias="schema")
kind: Literal["adsorption", "surface", "defect"]
members: list[Member]
coefficients: dict[str, float]
normalisation: dict[str, float]
unit: Literal["eV", "J/m2"]
formula: str
spec: dict[str, Any]
reservoirs: dict[str, Reservoir] = Field(default_factory=dict)
assumptions: list[str]
references: list[str]
def _formula_of(symbols: list[str]) -> dict[str, int]:
return dict(sorted(Counter(symbols).items()))
def _structure(document: Document, identifier: str) -> PeriodicStructure:
found = {s.id: s for s in document.periodic_structures}
if identifier not in found:
raise ValueError(f"periodic structure {identifier!r} not found")
return found[identifier]
def _species(document: Document, identifier: str) -> Species:
found = {s.id: s for s in document.species}
if identifier not in found:
raise ValueError(f"species {identifier!r} not found")
species = found[identifier]
if not species.conformers:
raise ValueError(f"species {identifier!r} needs a 3D structure first (embed it)")
return species
def _method_label(spec: PeriodicJobSpec, multiplicity: int = 1) -> str:
"""The method a result QVF must state, for the evaluation check."""
if spec.model == "seccm":
return "msindo"
method = spec.method.lower()
return ("uhf" if multiplicity > 1 else "rhf") if method == "hf" else method
def _molecule_input(species: Species, spec: PeriodicJobSpec, stem: str) -> InputFile:
"""A gas-phase reference computed with the members' Hamiltonian."""
from vibesketcher import __version__
if spec.model not in ("pbc", "seccm"):
raise WriterRefusal(
"vibeqc-energy", [f"no matching gas-phase route for model {spec.model}"]
)
conformer = species.conformers[0]
rows = [
f" ({ATOMIC_NUMBER[s]:>3}, " + ", ".join(f"{x: .10f}" for x in xyz) + "),"
for s, xyz in zip(conformer.symbols, conformer.xyz_angstrom, strict=True)
]
if spec.model == "seccm":
call = [' method="msindo",']
else:
method = spec.method.lower()
open_shell = species.multiplicity > 1
if method == "hf":
call = [f' method="{"uhf" if open_shell else "rhf"}",']
else:
call = [
f' method="{"uks" if open_shell else "rks"}",',
f" functional={method!r},",
]
call.append(f" basis={spec.basis!r},")
lines = [
f'"""{species.id}: gas-phase reference for a solid-state energy, written by '
f'vibe-sketcher {__version__}."""',
"",
"import vibeqc as vq",
"",
f"ANGSTROM_TO_BOHR = {ANGSTROM_TO_BOHR!r}",
"ATOMS = [",
*rows,
"]",
"molecule = vq.Molecule(",
" [vq.Atom(z, [x * ANGSTROM_TO_BOHR for x in xyz]) for z, *xyz in ATOMS],",
f" charge={species.charge},",
f" multiplicity={species.multiplicity},",
")",
"result = vq.run_job(",
" molecule,",
*call,
f" output={stem!r},",
" output_qvf=True,",
")",
'print(f"E = {result.energy:.10f} Ha")',
]
return InputFile(f"{stem}.py", "\n".join(lines) + "\n", "text/x-python")
def _member_input(structure: PeriodicStructure, spec: PeriodicJobSpec, role: str) -> InputFile:
bundle = write_periodic_input(structure, spec)
file = bundle.files[0]
stem = f"{role}-{file.name}"
content = file.content.replace(f"output={file.name[:-3]!r}", f"output={stem[:-3]!r}")
return InputFile(stem, content, file.media_type)
def _area(structure: PeriodicStructure) -> float:
rows = crystal.lattice_rows(structure)
return float(np.linalg.norm(np.cross(rows[0], rows[1])))
[docs]
def plan_adsorption(
document: Document, complex_id: str, slab_id: str, gas: str, spec: PeriodicJobSpec
) -> tuple[SolidEnergyPlan, InputBundle]:
complex_, slab = _structure(document, complex_id), _structure(document, slab_id)
builder = crystal._builder(complex_)
if builder.get("kind") != "adsorbate" or builder.get("slab") != slab_id:
raise ValueError(f"{complex_id} was not built by placing adsorbates on {slab_id}")
if not np.allclose(crystal.lattice_rows(complex_)[:2], crystal.lattice_rows(slab)[:2]):
raise ValueError("the complex and the clean slab have different in-plane cells")
molecule = _species(document, gas)
gas_formula = Counter(molecule.conformers[0].symbols)
excess = Counter(s.symbol for s in complex_.sites)
excess.subtract(Counter(s.symbol for s in slab.sites))
excess = Counter({k: v for k, v in excess.items() if v})
counts = {excess[k] / v for k, v in gas_formula.items()}
if (
set(excess) != set(gas_formula)
or len(counts) != 1
or not float(next(iter(counts))).is_integer()
):
raise ValueError(
f"{complex_id} minus {slab_id} is {dict(excess)}, not a whole number of {gas}"
)
n = int(next(iter(counts)))
files = [
_member_input(complex_, spec, "complex"),
_member_input(slab, spec, "slab"),
_molecule_input(molecule, spec, f"gas-{gas}"),
]
members = [
Member(
role="complex",
structure=complex_id,
formula=_formula_of([s.symbol for s in complex_.sites]),
input=files[0].name,
method=_method_label(spec),
),
Member(
role="slab",
structure=slab_id,
formula=_formula_of([s.symbol for s in slab.sites]),
input=files[1].name,
method=_method_label(spec),
),
Member(
role="gas",
species=gas,
formula=dict(sorted(gas_formula.items())),
input=files[2].name,
method=_method_label(spec, molecule.multiplicity),
),
]
plan = SolidEnergyPlan(
kind="adsorption",
members=members,
coefficients={"complex": 1 / n, "slab": -1 / n, "gas": -1.0},
normalisation={"adsorbates": n, "area_angstrom2": _area(slab)},
unit="eV",
formula=f"E_ads = [E(complex) - E(slab) - {n} E({gas})] / {n}",
spec=spec.model_dump(mode="json"),
assumptions=[
"the gas-phase molecule is computed as an isolated molecule with the same method "
"and basis; no counterpoise correction",
"geometries are used as built; relax them first for adsorption energies",
f"coverage {n} adsorbate(s) per {_area(slab):.3f} A^2 surface cell",
],
references=[],
)
return plan, InputBundle("vibeqc", complex_id, tuple(files))
[docs]
def plan_surface(
document: Document,
slab_id: str,
bulk_id: str,
spec: PeriodicJobSpec,
bulk_spec: PeriodicJobSpec | None = None,
) -> tuple[SolidEnergyPlan, InputBundle]:
slab, bulk = _structure(document, slab_id), _structure(document, bulk_id)
builder = crystal._builder(slab)
if builder.get("kind") != "slab":
raise ValueError(f"{slab_id} is not a clean slab made by the slab builder")
warnings = builder.get("warnings") or []
if warnings:
raise ValueError(
"a surface energy from one slab needs a stoichiometric slab with identical faces: "
+ "; ".join(str(w) for w in warnings)
)
if sum(bulk.pbc) != 3:
raise ValueError("the bulk reference must be three-dimensional")
slab_counts = Counter(s.symbol for s in slab.sites)
bulk_counts = Counter(s.symbol for s in bulk.sites)
ratios = {slab_counts[k] / v for k, v in bulk_counts.items()}
if set(slab_counts) != set(bulk_counts) or len(ratios) != 1:
raise ValueError(f"{slab_id} is not a whole multiple of the bulk cell {bulk_id}")
multiple = float(next(iter(ratios)))
area = _area(slab)
files = [_member_input(slab, spec, "slab"), _member_input(bulk, bulk_spec or spec, "bulk")]
members = [
Member(
role="slab",
structure=slab_id,
formula=dict(sorted(slab_counts.items())),
input=files[0].name,
method=_method_label(spec),
),
Member(
role="bulk",
structure=bulk_id,
formula=dict(sorted(bulk_counts.items())),
input=files[1].name,
method=_method_label(bulk_spec or spec),
),
]
scale = 1 / (2 * area)
plan = SolidEnergyPlan(
kind="surface",
members=members,
coefficients={"slab": scale, "bulk": -multiple * scale},
normalisation={"area_angstrom2": area, "bulk_cells": multiple, "faces": 2},
unit="J/m2",
formula=f"gamma = [E(slab) - {multiple:g} E(bulk cell)] / (2 x {area:.4f} A^2)",
spec=spec.model_dump(mode="json"),
assumptions=[
"both faces are identical and unrelaxed or relaxed alike",
"the bulk energy per cell is computed with the same method and basis; converge the "
"bulk k-mesh separately",
"slab thickness and vacuum are not converged by the plan",
],
references=["J. C. Boettger, Phys. Rev. B 49, 16798 (1994)"],
)
return plan, InputBundle("vibeqc", slab_id, tuple(files))
[docs]
def plan_defect(
document: Document,
defective_id: str,
pristine_id: str,
reservoirs: dict[str, Reservoir],
spec: PeriodicJobSpec,
) -> tuple[SolidEnergyPlan, InputBundle]:
defective, pristine = _structure(document, defective_id), _structure(document, pristine_id)
builder = crystal._builder(defective)
if builder.get("pristine") != pristine_id:
raise ValueError(f"{defective_id} was not made from {pristine_id} by the defect builder")
if defective.charge not in (None, 0) or (spec.charge not in (None, 0)):
raise ValueError(
"charged defects need a Fermi level, potential alignment and a finite-size "
"correction (Freysoldt et al. 2009); only neutral defects are evaluated"
)
change = builder.get("stoichiometry_change")
if not isinstance(change, dict) or not change:
raise ValueError("the defect records no change of composition")
missing = sorted(set(change) - set(reservoirs))
if missing:
raise ValueError(f"state a reservoir (chemical potential) for {', '.join(missing)}")
files = [_member_input(defective, spec, "defective"), _member_input(pristine, spec, "pristine")]
members = [
Member(
role="defective",
structure=defective_id,
formula=_formula_of([s.symbol for s in defective.sites]),
input=files[0].name,
method=_method_label(spec),
),
Member(
role="pristine",
structure=pristine_id,
formula=_formula_of([s.symbol for s in pristine.sites]),
input=files[1].name,
method=_method_label(spec),
),
]
coefficients = {"defective": 1.0, "pristine": -1.0}
terms = []
for element, delta in sorted(change.items()):
reservoir = reservoirs[element]
if reservoir.species:
molecule = _species(document, reservoir.species)
atoms = Counter(molecule.conformers[0].symbols)[element]
if atoms != reservoir.atoms:
raise ValueError(
f"{reservoir.species} has {atoms} {element} atoms, not {reservoir.atoms}"
)
role = f"mu-{element}"
file = _molecule_input(molecule, spec, f"{role}-{reservoir.species}")
files.append(file)
members.append(
Member(
role=role,
species=reservoir.species,
formula=dict(sorted(Counter(molecule.conformers[0].symbols).items())),
input=file.name,
method=_method_label(spec, molecule.multiplicity),
)
)
coefficients[role] = -float(delta) / reservoir.atoms
terms.append(f"{-float(delta):+g} E({reservoir.species})/{reservoir.atoms}")
elif reservoir.value_ev is not None and reservoir.source:
terms.append(f"{-float(delta):+g} x {reservoir.value_ev:g} eV ({reservoir.source})")
else:
raise ValueError(f"the {element} reservoir needs value_ev with a source, or a species")
constant = sum(
-float(delta) * (reservoirs[el].value_ev or 0.0)
for el, delta in change.items()
if not reservoirs[el].species
)
plan = SolidEnergyPlan(
kind="defect",
members=members,
coefficients=coefficients,
normalisation={"constant_ev": constant},
unit="eV",
formula="E_f = E(defective) - E(pristine) " + " ".join(terms),
spec=spec.model_dump(mode="json"),
reservoirs=reservoirs,
assumptions=[
"neutral defect; no charge corrections apply",
"defect-defect interaction through the periodic images is not removed; converge the "
"supercell size",
],
references=[
"C. Freysoldt, J. Neugebauer and C. G. Van de Walle, Phys. Rev. Lett. 102, "
"016402 (2009) (why charged defects are refused)"
],
)
return plan, InputBundle("vibeqc", defective_id, tuple(files))
def _read_result(data: bytes) -> dict[str, Any]:
"""Energy, method and composition from a vibe-qc result QVF."""
from vibesketcher.limits import check_archive
check_archive(data)
with zipfile.ZipFile(io.BytesIO(data)) as archive:
names = archive.namelist()
manifest = json.loads(archive.read(next(n for n in names if n.endswith("manifest.json"))))
provenance = manifest.get("provenance", {})
structures = [s for s in manifest.get("sections", []) if s.get("kind") == "structure"]
symbols: list[str] = []
if structures:
section = structures[0]
member = section.get("members", {}).get("structure", {})
path = member.get("path") if isinstance(member, dict) else None
prefix = next(
(n[: -len("manifest.json")] for n in names if n.endswith("manifest.json")), ""
)
if path:
payload = json.loads(archive.read(prefix + path))
symbols = [str(a.get("symbol")) for a in payload.get("atoms", [])]
energy = provenance.get("scf_energy") or {}
if not isinstance(energy, dict) or str(energy.get("units", "")).lower() not in (
"hartree",
"eh",
"ha",
):
raise ValueError("the result has no SCF energy in hartree")
status = str(provenance.get("run_status", "completed"))
if (
status not in ("completed", "success", "converged")
or provenance.get("scf_converged") is False
):
raise ValueError(f"the run is {status}, not a converged result")
runtime = (provenance.get("runtime") or {}).get("vibeqc") or {}
return {
"energy_hartree": float(energy["value"]),
"method": str(provenance.get("method", "")),
"formula": _formula_of(symbols) if symbols else None,
"vibeqc": {k: runtime.get(k) for k in ("version", "git_sha") if runtime.get(k)},
}
_METHOD_ALIASES = {"msindo": {"msindo", "ccm", "seccm"}, "rhf": {"rhf", "hf"}, "uhf": {"uhf"}}
[docs]
def evaluate(plan: SolidEnergyPlan, results: dict[str, bytes]) -> dict[str, Any]:
"""Combine member results into the planned energy, refusing mismatched results."""
missing = sorted({m.role for m in plan.members} - set(results))
if missing:
raise ValueError(f"results missing for {', '.join(missing)}")
energies: dict[str, float] = {}
runtimes: dict[str, dict[str, str]] = {}
for member in plan.members:
result = _read_result(results[member.role])
if result["formula"] is not None and result["formula"] != member.formula:
raise ValueError(
f"the {member.role} result is {result['formula']}, not the planned {member.formula}"
)
stated = re.sub(r"[^a-z0-9]", "", result["method"].lower())
allowed = _METHOD_ALIASES.get(member.method, {member.method, "rks", "uks"})
if stated and stated not in {re.sub(r"[^a-z0-9]", "", a) for a in allowed}:
raise ValueError(
f"the {member.role} result was computed with {result['method']}, "
f"not {member.method}"
)
energies[member.role] = result["energy_hartree"]
runtimes[member.role] = result["vibeqc"]
total_ha = sum(plan.coefficients[role] * energies[role] for role in plan.coefficients)
value_ev = total_ha * HARTREE_EV + plan.normalisation.get("constant_ev", 0.0)
report: dict[str, Any] = {
"kind": plan.kind,
"formula": plan.formula,
"member_energies_hartree": energies,
"vibeqc": runtimes,
"assumptions": plan.assumptions,
"references": plan.references,
}
if plan.unit == "J/m2":
report["value"] = value_ev * EV_PER_A2_TO_J_PER_M2
report["unit"] = "J/m2"
report["value_ev_per_angstrom2"] = value_ev
else:
report["value"] = value_ev
report["unit"] = "eV"
report["value_kj_per_mol"] = value_ev * EV_KJMOL
if not math.isfinite(report["value"]):
raise ValueError("the energy is not finite")
return report