Skills › Scientific Agent Skills › pymatgen › phase_diagram_generator.py
phase_diagram_generator.py
416 lines · skills/pymatgen/scripts/phase_diagram_generator.py at 5f2bc03
1#!/usr/bin/env python3
2"""Build a local phase diagram from a strict, provenance-bearing JSON dataset."""
4from __future__ import annotations
6import argparse
7import math
8import tempfile
9from pathlib import Path
10from typing import Any
12from _common import (
13 ABSOLUTE_MAX_OUTPUT_BYTES,
14 CliError,
15 DEFAULT_MAX_INPUT_BYTES,
16 DEFAULT_MAX_OUTPUT_BYTES,
17 atomic_link_from_temp,
18 checked_output_file,
19 emit_json,
20 finite_float,
21 load_strict_json,
22 positive_int,
23 write_json_new,
24)
27TOP_LEVEL_KEYS = {
28 "schema_version",
29 "energy_unit",
30 "energy_basis",
31 "provenance",
32 "entries",
33}
34PROVENANCE_KEYS = {
35 "source",
36 "method",
37 "retrieved_at",
38 "database_version",
39 "license",
40 "citation",
41 "notes",
42}
43ENTRY_KEYS = {"entry_id", "composition", "energy_eV", "provenance"}
46def validate_provenance(value: Any, label: str) -> dict[str, str]:
47 """Validate a small string-only provenance object."""
48 if not isinstance(value, dict):
49 raise CliError(f"{label} must be an object")
50 unknown = set(value) - PROVENANCE_KEYS
51 if unknown:
52 raise CliError(f"{label} has unknown keys: {sorted(unknown)}")
53 if not isinstance(value.get("source"), str) or not value["source"].strip():
54 raise CliError(f"{label}.source must be a non-empty string")
55 result: dict[str, str] = {}
56 for key, item in value.items():
57 if not isinstance(item, str) or len(item) > 2000:
58 raise CliError(f"{label}.{key} must be a string of at most 2000 chars")
59 result[key] = item
60 return result
63def entries_from_payload(
64 payload: Any,
65 *,
66 max_entries: int,
67) -> tuple[list[Any], dict[str, Any]]:
68 """Validate the strict schema and construct plain ComputedEntry objects."""
69 if not isinstance(payload, dict):
70 raise CliError("top-level JSON value must be an object")
71 missing = TOP_LEVEL_KEYS - set(payload)
72 unknown = set(payload) - TOP_LEVEL_KEYS
73 if missing or unknown:
74 raise CliError(
75 f"top-level keys mismatch; missing={sorted(missing)}, "
76 f"unknown={sorted(unknown)}"
77 )
78 if payload["schema_version"] != "1.0":
79 raise CliError("schema_version must be exactly '1.0'")
80 if payload["energy_unit"] != "eV":
81 raise CliError("energy_unit must be exactly 'eV'")
82 if payload["energy_basis"] != "total_per_entry":
83 raise CliError("energy_basis must be exactly 'total_per_entry'")
84 dataset_provenance = validate_provenance(
85 payload["provenance"], "provenance"
86 )
87 rows = payload["entries"]
88 if not isinstance(rows, list) or not rows:
89 raise CliError("entries must be a non-empty array")
90 if len(rows) > max_entries:
91 raise CliError(f"entries exceeds the {max_entries}-entry limit")
93 from pymatgen.core import Composition
94 from pymatgen.entries.computed_entries import ComputedEntry
96 entries: list[Any] = []
97 seen_ids: set[str] = set()
98 provenance_by_id: dict[str, dict[str, str]] = {}
99 for index, row in enumerate(rows):
100 label = f"entries[{index}]"
101 if not isinstance(row, dict) or set(row) != ENTRY_KEYS:
102 actual = sorted(row) if isinstance(row, dict) else type(row).__name__
103 raise CliError(f"{label} must contain exactly {sorted(ENTRY_KEYS)}; got {actual}")
104 entry_id = row["entry_id"]
105 if (
106 not isinstance(entry_id, str)
107 or not entry_id.strip()
108 or len(entry_id) > 200
109 ):
110 raise CliError(f"{label}.entry_id must be a short non-empty string")
111 if entry_id in seen_ids:
112 raise CliError(f"duplicate entry_id: {entry_id!r}")
113 seen_ids.add(entry_id)
114 formula = row["composition"]
115 if not isinstance(formula, str) or len(formula) > 500:
116 raise CliError(f"{label}.composition must be a formula string")
117 try:
118 composition = Composition(formula, strict=True)
119 except (TypeError, ValueError) as exc:
120 raise CliError(f"{label}.composition is invalid: {exc}") from exc
121 if composition.num_atoms <= 0:
122 raise CliError(f"{label}.composition must contain positive amounts")
123 energy = finite_float(row["energy_eV"], f"{label}.energy_eV")
124 entry_provenance = validate_provenance(
125 row["provenance"], f"{label}.provenance"
126 )
127 provenance_by_id[entry_id] = entry_provenance
128 entries.append(
129 ComputedEntry(
130 composition,
131 energy,
132 entry_id=entry_id,
133 data={"provenance": entry_provenance},
134 )
135 )
136 return entries, {
137 "dataset": dataset_provenance,
138 "entries": provenance_by_id,
139 }
142def phase_report(
143 entries: list[Any],
144 provenance: dict[str, Any],
145 *,
146 analyze: list[str],
147 max_report_entries: int,
148) -> tuple[dict[str, Any], Any]:
149 """Construct a phase diagram and a bounded scientific report."""
150 from pymatgen.analysis.phase_diagram import PhaseDiagram
151 from pymatgen.core import Composition
153 diagram = PhaseDiagram(entries)
154 stable = set(diagram.stable_entries)
155 ordered_entries = sorted(
156 entries,
157 key=lambda item: (
158 item.composition.reduced_formula,
159 float(item.energy_per_atom),
160 str(item.entry_id),
161 ),
162 )
163 rows: list[dict[str, Any]] = []
164 for entry in ordered_entries[:max_report_entries]:
165 rows.append(
166 {
167 "entry_id": str(entry.entry_id),
168 "formula": entry.composition.reduced_formula,
169 "energy_eV_total": float(entry.energy),
170 "energy_eV_per_atom": float(entry.energy_per_atom),
171 "formation_energy_eV_per_atom": float(
172 diagram.get_form_energy_per_atom(entry)
173 ),
174 "energy_above_hull_eV_per_atom": float(
175 diagram.get_e_above_hull(entry)
176 ),
177 "on_computed_convex_hull": entry in stable,
178 }
179 )
181 analyses: list[dict[str, Any]] = []
182 for formula in analyze:
183 try:
184 composition = Composition(formula, strict=True)
185 decomposition = diagram.get_decomposition(composition)
186 matches = [
187 entry
188 for entry in entries
189 if entry.composition.fractional_composition
190 == composition.fractional_composition
191 ]
192 analyses.append(
193 {
194 "query": formula,
195 "reduced_formula": composition.reduced_formula,
196 "matching_entries": [
197 {
198 "entry_id": str(entry.entry_id),
199 "energy_above_hull_eV_per_atom": float(
200 diagram.get_e_above_hull(entry)
201 ),
202 }
203 for entry in sorted(
204 matches,
205 key=lambda item: (
206 float(item.energy_per_atom),
207 str(item.entry_id),
208 ),
209 )
210 ],
211 "computed_hull_decomposition": [
212 {
213 "entry_id": str(entry.entry_id),
214 "formula": entry.composition.reduced_formula,
215 "fraction": float(fraction),
216 }
217 for entry, fraction in sorted(
218 decomposition.items(),
219 key=lambda item: str(item[0].entry_id),
220 )
221 ],
222 }
223 )
224 except (TypeError, ValueError) as exc:
225 analyses.append(
226 {
227 "query": formula,
228 "error": f"{type(exc).__name__}: {exc}"[:500],
229 }
230 )
232 report = {
233 "ok": True,
234 "analysis": "local_computed_phase_diagram",
235 "units": {
236 "input_energy": "eV total per entry",
237 "reported_normalized_energy": "eV/atom",
238 },
239 "chemical_system": "-".join(str(element) for element in diagram.elements),
240 "elements": [str(element) for element in diagram.elements],
241 "entry_count": len(entries),
242 "stable_entry_count": len(stable),
243 "entries": rows,
244 "entries_omitted": max(0, len(entries) - max_report_entries),
245 "composition_analyses": analyses,
246 "provenance": provenance["dataset"],
247 "interpretation_limits": [
248 "The hull is conditional on this exact entry set and energy model.",
249 "Energies from incompatible methods or correction schemes must not be mixed.",
250 "Computed stability is not experimental truth or a synthesis guarantee.",
251 "Finite-temperature, pressure, kinetic, disorder, and uncertainty effects are absent unless encoded upstream.",
252 ],
253 "experimental_validity_established": False,
254 }
255 return report, diagram
258def write_plot_new(
259 diagram: Any,
260 output_path: Path,
261 *,
262 show_unstable: float,
263 max_bytes: int,
264) -> None:
265 """Render to a sibling temporary file, then link without overwriting."""
266 suffix = output_path.suffix.casefold()
267 if suffix not in {".png", ".pdf", ".svg"}:
268 raise CliError("plot output suffix must be .png, .pdf, or .svg")
269 from pymatgen.analysis.phase_diagram import PDPlotter
271 temporary_path: Path | None = None
272 try:
273 with tempfile.NamedTemporaryFile(
274 prefix=".pymatgen-phase-",
275 suffix=suffix,
276 dir=output_path.parent,
277 delete=False,
278 ) as handle:
279 temporary_path = Path(handle.name)
280 # Matplotlib exports locally without Plotly's optional Kaleido/browser stack.
281 plotter = PDPlotter(diagram, show_unstable=show_unstable, backend="matplotlib")
282 plotter.write_image(
283 str(temporary_path),
284 image_format=suffix.removeprefix("."),
285 )
286 if temporary_path.stat().st_size > max_bytes:
287 raise CliError(f"plot exceeds the {max_bytes}-byte limit")
288 atomic_link_from_temp(temporary_path, output_path)
289 finally:
290 if temporary_path is not None:
291 temporary_path.unlink(missing_ok=True)
294def build_parser() -> argparse.ArgumentParser:
295 parser = argparse.ArgumentParser(
296 description=(
297 "Build an offline phase diagram from strict JSON total energies. "
298 "This script never queries Materials Project."
299 )
300 )
301 parser.add_argument("entries_json", help="Strict local entries dataset")
302 parser.add_argument(
303 "--analyze",
304 action="append",
305 default=[],
306 help="Formula to decompose or match (repeatable, maximum 20)",
307 )
308 parser.add_argument("--output", help="New JSON report; default is stdout")
309 parser.add_argument("--plot", help="New .png, .pdf, or .svg plot")
310 parser.add_argument(
311 "--show-unstable",
312 type=float,
313 default=0.2,
314 help="Plot unstable entries up to this eV/atom (default: 0.2)",
315 )
316 parser.add_argument(
317 "--max-input-bytes",
318 type=positive_int,
319 default=DEFAULT_MAX_INPUT_BYTES,
320 )
321 parser.add_argument("--max-entries", type=positive_int, default=5000)
322 parser.add_argument("--max-report-entries", type=positive_int, default=200)
323 parser.add_argument(
324 "--max-output-bytes",
325 type=positive_int,
326 default=DEFAULT_MAX_OUTPUT_BYTES,
327 )
328 return parser
331def main() -> int:
332 args = build_parser().parse_args()
333 try:
334 if len(args.analyze) > 20:
335 raise CliError("--analyze may be repeated at most 20 times")
336 if not math.isfinite(args.show_unstable) or args.show_unstable < 0:
337 raise CliError("--show-unstable must be finite and non-negative")
338 if args.max_entries > 10_000:
339 raise CliError("--max-entries may not exceed 10000")
340 if args.max_report_entries > 1_000:
341 raise CliError("--max-report-entries may not exceed 1000")
342 if args.max_output_bytes > ABSOLUTE_MAX_OUTPUT_BYTES:
343 raise CliError(
344 f"--max-output-bytes may not exceed {ABSOLUTE_MAX_OUTPUT_BYTES}"
345 )
346 payload, input_path = load_strict_json(
347 args.entries_json,
348 max_bytes=args.max_input_bytes,
349 )
350 entries, provenance = entries_from_payload(
351 payload,
352 max_entries=args.max_entries,
353 )
354 report, diagram = phase_report(
355 entries,
356 provenance,
357 analyze=args.analyze,
358 max_report_entries=args.max_report_entries,
359 )
360 report["input"] = {
361 "name": input_path.name,
362 "bytes": input_path.stat().st_size,
363 "network_accessed": False,
364 }
365 created: list[str] = []
366 if args.output:
367 output_path = checked_output_file(
368 args.output,
369 input_paths=(input_path,),
370 )
371 write_json_new(
372 output_path,
373 report,
374 max_bytes=args.max_output_bytes,
375 )
376 created.append(output_path.name)
377 if args.plot:
378 plot_path = checked_output_file(
379 args.plot,
380 input_paths=(input_path,),
381 )
382 if len(diagram.elements) > 4:
383 raise CliError("plotting is limited to at most four elements")
384 write_plot_new(
385 diagram,
386 plot_path,
387 show_unstable=args.show_unstable,
388 max_bytes=args.max_output_bytes,
389 )
390 created.append(plot_path.name)
391 if created:
392 emit_json(
393 {
394 "ok": True,
395 "created": created,
396 "overwrote_existing": False,
397 "network_accessed": False,
398 "entry_count": len(entries),
399 }
400 )
401 else:
402 emit_json(report)
403 return 0
404 except (CliError, ImportError, OSError, RuntimeError, TypeError, ValueError) as exc:
405 emit_json(
406 {
407 "ok": False,
408 "error": f"{type(exc).__name__}: {exc}"[:1000],
409 "network_accessed": False,
410 }
411 )
412 return 2
415if __name__ == "__main__":
416 raise SystemExit(main())