Skills › Scientific Agent Skills › pymatgen › structure_analyzer.py
structure_analyzer.py
294 lines · skills/pymatgen/scripts/structure_analyzer.py at 5f2bc03
1#!/usr/bin/env python3
2"""Produce a bounded JSON analysis of one local periodic structure."""
4from __future__ import annotations
6import argparse
7import math
8import warnings
9from typing import Any
11from _common import (
12 ABSOLUTE_MAX_PAIRWISE_SITES,
13 CliError,
14 DEFAULT_MAX_INPUT_BYTES,
15 DEFAULT_MAX_OUTPUT_BYTES,
16 DEFAULT_MAX_SITES,
17 checked_output_file,
18 emit_json,
19 load_structure,
20 positive_int,
21 structure_oxidation_summary,
22 write_json_new,
23)
26def minimum_periodic_distance(structure: Any, max_sites: int) -> float | None:
27 """Return the minimum non-self periodic distance under an explicit bound."""
28 if len(structure) < 2 or len(structure) > max_sites:
29 return None
30 minimum = math.inf
31 for first in range(len(structure)):
32 for second in range(first + 1, len(structure)):
33 minimum = min(minimum, float(structure.get_distance(first, second)))
34 return minimum if math.isfinite(minimum) else None
37def site_records(structure: Any, limit: int) -> list[dict[str, Any]]:
38 """Return a bounded, explicit fractional-coordinate site table."""
39 records: list[dict[str, Any]] = []
40 for index, site in enumerate(structure[:limit]):
41 records.append(
42 {
43 "index": index,
44 "species": {
45 str(specie): float(occupancy)
46 for specie, occupancy in site.species.items()
47 },
48 "fractional_coordinates": [
49 float(value) for value in site.frac_coords
50 ],
51 "label": site.label,
52 }
53 )
54 return records
57def symmetry_report(
58 structure: Any,
59 *,
60 symprec: float,
61 angle_tolerance: float,
62) -> dict[str, Any]:
63 """Run one explicitly parameterized spglib-backed symmetry analysis."""
64 from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
66 analyzer = SpacegroupAnalyzer(
67 structure,
68 symprec=symprec,
69 angle_tolerance=angle_tolerance,
70 )
71 symmetrized = analyzer.get_symmetrized_structure()
72 return {
73 "backend": "spglib through pymatgen",
74 "symprec_angstrom": symprec,
75 "angle_tolerance_degrees": angle_tolerance,
76 "space_group_symbol": analyzer.get_space_group_symbol(),
77 "space_group_number": analyzer.get_space_group_number(),
78 "crystal_system": str(analyzer.get_crystal_system()),
79 "point_group_symbol": analyzer.get_point_group_symbol(),
80 "symmetry_operations": len(analyzer.get_symmetry_operations()),
81 "equivalent_site_groups": len(symmetrized.equivalent_indices),
82 "wyckoff_symbols_by_group": list(symmetrized.wyckoff_symbols),
83 "tolerance_sensitivity_assessed": False,
84 }
87def neighbor_report(
88 structure: Any,
89 *,
90 site_limit: int,
91 neighbor_limit: int,
92) -> dict[str, Any]:
93 """Run CrystalNN for a bounded site prefix and truncate each neighbor list."""
94 if not structure.is_ordered:
95 raise CliError("CrystalNN report requires an ordered structure")
96 from pymatgen.analysis.local_env import CrystalNN
98 strategy = CrystalNN()
99 rows: list[dict[str, Any]] = []
100 warning_messages: list[str] = []
101 with warnings.catch_warnings(record=True) as caught:
102 warnings.simplefilter("always")
103 for index in range(min(len(structure), site_limit)):
104 try:
105 neighbors = strategy.get_nn_info(structure, index)
106 rows.append(
107 {
108 "site_index": index,
109 "species": structure[index].species_string,
110 "coordination_number": len(neighbors),
111 "neighbors": [
112 {
113 "site_index": int(item["site_index"]),
114 "image": [int(value) for value in item["image"]],
115 "weight": float(item["weight"]),
116 }
117 for item in neighbors[:neighbor_limit]
118 ],
119 "neighbors_omitted": max(
120 0, len(neighbors) - neighbor_limit
121 ),
122 }
123 )
124 except (RuntimeError, TypeError, ValueError) as exc:
125 rows.append(
126 {
127 "site_index": index,
128 "error": f"{type(exc).__name__}: {exc}"[:500],
129 }
130 )
131 warning_messages = [
132 f"{item.category.__name__}: {item.message}" for item in caught[:20]
133 ]
134 return {
135 "method": "CrystalNN",
136 "sites": rows,
137 "sites_omitted": max(0, len(structure) - site_limit),
138 "warnings": warning_messages,
139 "coordination_is_model_dependent": True,
140 }
143def analyze_structure(structure: Any, args: argparse.Namespace) -> dict[str, Any]:
144 """Build a bounded analysis payload."""
145 lattice = structure.lattice
146 report: dict[str, Any] = {
147 "ok": True,
148 "analysis": "periodic_structure",
149 "units": {
150 "length": "angstrom",
151 "angle": "degree",
152 "volume": "angstrom^3",
153 "density": "g/cm^3",
154 "mass": "amu per composition represented",
155 },
156 "periodicity": {
157 "periodic_boundary_conditions": [bool(value) for value in lattice.pbc],
158 "lattice_required": True,
159 },
160 "composition": {
161 "formula": structure.composition.formula,
162 "reduced_formula": structure.composition.reduced_formula,
163 "hill_formula": structure.composition.hill_formula,
164 "chemical_system": structure.composition.chemical_system,
165 "mass_amu": float(structure.composition.weight),
166 "charge": float(structure.charge),
167 "ordered": bool(structure.is_ordered),
168 "oxidation_states": structure_oxidation_summary(structure),
169 },
170 "lattice": {
171 "matrix_rows_angstrom": [
172 [float(value) for value in row] for row in lattice.matrix
173 ],
174 "abc_angstrom": [float(value) for value in lattice.abc],
175 "angles_degrees": [float(value) for value in lattice.angles],
176 "volume_angstrom_cubed": float(structure.volume),
177 "density_g_cm3": float(structure.density),
178 },
179 "sites": {
180 "count": len(structure),
181 "coordinate_mode": "fractional",
182 "records": site_records(structure, args.max_site_records),
183 "records_omitted": max(0, len(structure) - args.max_site_records),
184 },
185 "minimum_periodic_distance_angstrom": minimum_periodic_distance(
186 structure, args.max_distance_sites
187 ),
188 "minimum_distance_omitted_above_sites": args.max_distance_sites,
189 "minimum_distance_includes_self_images": False,
190 "scientific_validity_established": False,
191 }
192 if args.symmetry:
193 report["symmetry"] = symmetry_report(
194 structure,
195 symprec=args.symprec,
196 angle_tolerance=args.angle_tolerance,
197 )
198 if args.neighbors:
199 report["neighbors"] = neighbor_report(
200 structure,
201 site_limit=args.max_neighbor_sites,
202 neighbor_limit=args.max_neighbors_per_site,
203 )
204 return report
207def build_parser() -> argparse.ArgumentParser:
208 parser = argparse.ArgumentParser(
209 description=(
210 "Analyze one bounded local periodic structure. JSON output is "
211 "truncated by explicit site/neighbor limits."
212 )
213 )
214 parser.add_argument("structure_file", help="Existing local structure file")
215 parser.add_argument("--structure-index", type=int, default=0)
216 parser.add_argument("--symmetry", action="store_true")
217 parser.add_argument("--symprec", type=float, default=0.01)
218 parser.add_argument("--angle-tolerance", type=float, default=5.0)
219 parser.add_argument("--neighbors", action="store_true")
220 parser.add_argument("--output", help="New JSON output; default is stdout")
221 parser.add_argument(
222 "--max-input-bytes",
223 type=positive_int,
224 default=DEFAULT_MAX_INPUT_BYTES,
225 )
226 parser.add_argument("--max-sites", type=positive_int, default=DEFAULT_MAX_SITES)
227 parser.add_argument("--max-site-records", type=positive_int, default=100)
228 parser.add_argument("--max-distance-sites", type=positive_int, default=500)
229 parser.add_argument("--max-neighbor-sites", type=positive_int, default=100)
230 parser.add_argument(
231 "--max-neighbors-per-site", type=positive_int, default=24
232 )
233 parser.add_argument(
234 "--max-output-bytes",
235 type=positive_int,
236 default=DEFAULT_MAX_OUTPUT_BYTES,
237 )
238 return parser
241def main() -> int:
242 args = build_parser().parse_args()
243 try:
244 if args.structure_index < 0:
245 raise CliError("--structure-index must be non-negative")
246 if not math.isfinite(args.symprec) or args.symprec <= 0:
247 raise CliError("--symprec must be finite and positive")
248 if not math.isfinite(args.angle_tolerance) or args.angle_tolerance < 0:
249 raise CliError("--angle-tolerance must be finite and non-negative")
250 if args.max_distance_sites > ABSOLUTE_MAX_PAIRWISE_SITES:
251 raise CliError(
252 f"--max-distance-sites may not exceed "
253 f"{ABSOLUTE_MAX_PAIRWISE_SITES}"
254 )
255 if args.max_site_records > 10_000:
256 raise CliError("--max-site-records may not exceed 10000")
257 if args.max_neighbor_sites > 10_000:
258 raise CliError("--max-neighbor-sites may not exceed 10000")
259 if args.max_neighbors_per_site > 1_000:
260 raise CliError("--max-neighbors-per-site may not exceed 1000")
261 structure, input_path, parse_report = load_structure(
262 args.structure_file,
263 structure_index=args.structure_index,
264 max_bytes=args.max_input_bytes,
265 max_sites=args.max_sites,
266 )
267 report = analyze_structure(structure, args)
268 report["input"] = parse_report
269 if args.output:
270 output = checked_output_file(args.output, input_paths=(input_path,))
271 write_json_new(
272 output,
273 report,
274 max_bytes=args.max_output_bytes,
275 )
276 emit_json(
277 {
278 "ok": True,
279 "output": output.name,
280 "output_bytes": output.stat().st_size,
281 "overwrote_existing": False,
282 "sites": len(structure),
283 }
284 )
285 else:
286 emit_json(report)
287 return 0
288 except (CliError, ImportError, OSError, RuntimeError, TypeError, ValueError) as exc:
289 emit_json({"ok": False, "error": f"{type(exc).__name__}: {exc}"[:1000]})
290 return 2
293if __name__ == "__main__":
294 raise SystemExit(main())