Update defect-analysis

quanshengwu · quanshengwu · commit 1f7d5fc7d376 · 2026-03-11T05:03:50.000+08:00
diff --git a/README.md b/README.md
@@ -4,6 +4,8 @@
 
 Standalone skill for defect-focused DFT result analysis.
 
+Supports VASP, QE, and ABINIT-style pristine/defect inputs.
+
 ## Install
 
 ```bash
@@ -16,6 +18,8 @@ npx skills add chatmaterials/defect-analysis -g -y
 python3 -m py_compile scripts/*.py
 npx skills add . --list
 python3 scripts/analyze_defect_formation.py fixtures/pristine fixtures/defect --species O --delta -1 --mu -4.0 --json
+python3 scripts/analyze_defect_formation.py fixtures/qe/pristine fixtures/qe/defect --species O --delta -1 --mu -4.0 --json
+python3 scripts/analyze_defect_formation.py fixtures/abinit/pristine fixtures/abinit/defect --species O --delta -1 --mu -4.0 --json
 python3 scripts/analyze_defect_structure.py fixtures/pristine/POSCAR fixtures/defect/POSCAR --json
 python3 scripts/export_defect_report.py fixtures/pristine fixtures/defect --species O --delta -1 --mu -4.0
 python3 scripts/run_regression.py
diff --git a/SKILL.md b/SKILL.md
@@ -1,6 +1,6 @@
 ---
 name: "defect-analysis"
-description: "Use when the task is to analyze defect-related DFT results, including neutral defect formation energies, stoichiometric changes, structural distortion around a defect, and compact markdown reports from pristine and defect calculations."
+description: "Use when the task is to analyze defect-related DFT results, including neutral defect formation energies, stoichiometric changes, structural distortion around a defect, and compact markdown reports from pristine and defect calculations. Supports VASP, QE, and ABINIT-style inputs."
 ---
 
 # Defect Analysis
@@ -14,6 +14,12 @@ Use this skill for defect-focused post-processing rather than generic workflow s
 - quantify structural or volume change caused by a defect
 - write a compact defect-analysis report from finished calculations
 
+Supported backends:
+
+- VASP-like directories with `OUTCAR` and `POSCAR`
+- QE-like directories with `.out` and structural input blocks
+- ABINIT-like directories with `.abo` and `.abi`
+
 ## Use the bundled helpers
 
 - `scripts/analyze_defect_formation.py`
diff --git a/fixtures/abinit/defect/Fe.psp8 b/fixtures/abinit/defect/Fe.psp8
@@ -0,0 +1 @@
+Pseudo placeholder
diff --git a/fixtures/abinit/defect/O.psp8 b/fixtures/abinit/defect/O.psp8
@@ -0,0 +1 @@
+Pseudo placeholder
diff --git a/fixtures/abinit/defect/defect.abi b/fixtures/abinit/defect/defect.abi
@@ -0,0 +1,19 @@
+ecut 20.0
+ngkpt 6 6 6
+nshiftk 1
+shiftk 0.0 0.0 0.0
+tolvrs 1.0d-10
+occopt 1
+pseudos "Fe.psp8, O.psp8"
+natom 2
+ntypat 2
+znucl 26 8
+typat 1 2
+acell 3*9.5431161
+rprim
+ 1.0 0.0 0.0
+ 0.0 1.0 0.0
+ 0.0 0.0 1.0
+xred
+ 0.0 0.0 0.0
+ 0.5 0.5 0.0
diff --git a/fixtures/abinit/defect/run.abo b/fixtures/abinit/defect/run.abo
@@ -0,0 +1,2 @@
+ etotal -0.77950703
+ Calculation completed.
diff --git a/fixtures/abinit/pristine/Fe.psp8 b/fixtures/abinit/pristine/Fe.psp8
@@ -0,0 +1 @@
+Pseudo placeholder
diff --git a/fixtures/abinit/pristine/O.psp8 b/fixtures/abinit/pristine/O.psp8
@@ -0,0 +1 @@
+Pseudo placeholder
diff --git a/fixtures/abinit/pristine/pristine.abi b/fixtures/abinit/pristine/pristine.abi
@@ -0,0 +1,20 @@
+ecut 20.0
+ngkpt 6 6 6
+nshiftk 1
+shiftk 0.0 0.0 0.0
+tolvrs 1.0d-10
+occopt 1
+pseudos "Fe.psp8, O.psp8"
+natom 3
+ntypat 2
+znucl 26 8
+typat 1 2 2
+acell 3*9.4486306
+rprim
+ 1.0 0.0 0.0
+ 0.0 1.0 0.0
+ 0.0 0.0 1.0
+xred
+ 0.0 0.0 0.0
+ 0.5 0.5 0.0
+ 0.5 0.0 0.5
diff --git a/fixtures/abinit/pristine/run.abo b/fixtures/abinit/pristine/run.abo
@@ -0,0 +1,2 @@
+ etotal -1.000000
+ Calculation completed.
diff --git a/fixtures/qe/defect/pseudo/Fe.UPF b/fixtures/qe/defect/pseudo/Fe.UPF
@@ -0,0 +1 @@
+<UPF>Fe</UPF>
diff --git a/fixtures/qe/defect/pseudo/O.UPF b/fixtures/qe/defect/pseudo/O.UPF
@@ -0,0 +1 @@
+<UPF>O</UPF>
diff --git a/fixtures/qe/defect/relax.in b/fixtures/qe/defect/relax.in
@@ -0,0 +1,32 @@
+&CONTROL
+  calculation = 'relax',
+  prefix = 'defect',
+  pseudo_dir = './pseudo',
+  outdir = './tmp',
+/
+&SYSTEM
+  ibrav = 0,
+  nat = 2,
+  ntyp = 2,
+  ecutwfc = 60.0,
+  ecutrho = 480.0,
+  occupations = 'fixed',
+/
+&ELECTRONS
+  conv_thr = 1.0d-8,
+/
+ATOMIC_SPECIES
+Fe 55.845 Fe.UPF
+O  15.999 O.UPF
+
+CELL_PARAMETERS angstrom
+5.05 0.00 0.00
+0.00 5.05 0.00
+0.00 0.00 5.05
+
+ATOMIC_POSITIONS crystal
+Fe 0.00 0.00 0.00
+O  0.50 0.50 0.00
+
+K_POINTS automatic
+6 6 6 0 0 0
diff --git a/fixtures/qe/defect/relax.out b/fixtures/qe/defect/relax.out
@@ -0,0 +1,5 @@
+ Program PWSCF v.7.2 starts
+ convergence has been achieved
+ !    total energy              =     -1.39647425 Ry
+ End of BFGS Geometry Optimization
+ JOB DONE.
diff --git a/fixtures/qe/pristine/pseudo/Fe.UPF b/fixtures/qe/pristine/pseudo/Fe.UPF
@@ -0,0 +1 @@
+<UPF>Fe</UPF>
diff --git a/fixtures/qe/pristine/pseudo/O.UPF b/fixtures/qe/pristine/pseudo/O.UPF
@@ -0,0 +1 @@
+<UPF>O</UPF>
diff --git a/fixtures/qe/pristine/relax.in b/fixtures/qe/pristine/relax.in
@@ -0,0 +1,33 @@
+&CONTROL
+  calculation = 'relax',
+  prefix = 'pristine',
+  pseudo_dir = './pseudo',
+  outdir = './tmp',
+/
+&SYSTEM
+  ibrav = 0,
+  nat = 3,
+  ntyp = 2,
+  ecutwfc = 60.0,
+  ecutrho = 480.0,
+  occupations = 'fixed',
+/
+&ELECTRONS
+  conv_thr = 1.0d-8,
+/
+ATOMIC_SPECIES
+Fe 55.845 Fe.UPF
+O  15.999 O.UPF
+
+CELL_PARAMETERS angstrom
+5.00 0.00 0.00
+0.00 5.00 0.00
+0.00 0.00 5.00
+
+ATOMIC_POSITIONS crystal
+Fe 0.00 0.00 0.00
+O  0.50 0.50 0.00
+O  0.50 0.00 0.50
+
+K_POINTS automatic
+6 6 6 0 0 0
diff --git a/fixtures/qe/pristine/relax.out b/fixtures/qe/pristine/relax.out
@@ -0,0 +1,5 @@
+ Program PWSCF v.7.2 starts
+ convergence has been achieved
+ !    total energy              =     -1.83746612 Ry
+ End of BFGS Geometry Optimization
+ JOB DONE.
diff --git a/scripts/analyze_defect_formation.py b/scripts/analyze_defect_formation.py
@@ -4,24 +4,19 @@
 
 import argparse
 import json
-import re
 from pathlib import Path
 
-
-def read_energy(path: Path) -> float:
-    outcar = path / "OUTCAR" if path.is_dir() else path
-    text = outcar.read_text(errors="ignore")
-    matches = re.findall(r"TOTEN\s*=\s*([\-0-9.Ee+]+)", text)
-    if not matches:
-        raise SystemExit(f"No TOTEN found in {outcar}")
-    return float(matches[-1])
+from defect_io import read_energy
 
 
 def analyze(pristine: Path, defect: Path, species: str, delta: int, mu: float) -> dict[str, object]:
-    e_bulk = read_energy(pristine)
-    e_def = read_energy(defect)
+    backend_pristine, e_bulk = read_energy(pristine)
+    backend_defect, e_def = read_energy(defect)
+    if backend_pristine != backend_defect:
+        raise SystemExit("Pristine and defect states must use the same backend for direct comparison")
     e_form = e_def - e_bulk - delta * mu
     return {
+        "backend": backend_pristine,
         "pristine": str(pristine),
         "defect": str(defect),
         "species": species,
diff --git a/scripts/analyze_defect_structure.py b/scripts/analyze_defect_structure.py
@@ -6,29 +6,14 @@
 import json
 from pathlib import Path
 
-
-def read_poscar(path: Path) -> dict[str, object]:
-    lines = path.read_text().splitlines()
-    scale = float(lines[1].split()[0])
-    lattice = [[float(x) * scale for x in lines[i].split()] for i in range(2, 5)]
-    species = lines[5].split()
-    counts = [int(x) for x in lines[6].split()]
-    natoms = sum(counts)
-    return {"lattice": lattice, "species": species, "counts": counts, "natoms": natoms}
-
-
-def volume(cell: list[list[float]]) -> float:
-    a, b, c = cell
-    return (
-        a[0] * (b[1] * c[2] - b[2] * c[1])
-        - a[1] * (b[0] * c[2] - b[2] * c[0])
-        + a[2] * (b[0] * c[1] - b[1] * c[0])
-    )
+from defect_io import read_structure, volume
 
 
 def analyze(pristine: Path, defect: Path) -> dict[str, object]:
-    a = read_poscar(pristine)
-    b = read_poscar(defect)
+    backend_a, a = read_structure(pristine)
+    backend_b, b = read_structure(defect)
+    if backend_a != backend_b:
+        raise SystemExit("Pristine and defect structures must use the same backend for direct comparison")
     vol_a = volume(a["lattice"])
     vol_b = volume(b["lattice"])
     species_delta = {}
@@ -38,6 +23,7 @@ def analyze(pristine: Path, defect: Path) -> dict[str, object]:
         count_b = b["counts"][b["species"].index(specie)] if specie in b["species"] else 0
         species_delta[specie] = count_b - count_a
     return {
+        "backend": backend_a,
         "pristine": str(pristine),
         "defect": str(defect),
         "natoms_pristine": a["natoms"],
diff --git a/scripts/defect_io.py b/scripts/defect_io.py
@@ -0,0 +1,127 @@
+#!/usr/bin/env python3
+
+from __future__ import annotations
+
+import re
+from pathlib import Path
+
+
+RY_TO_EV = 13.605693009
+HARTREE_TO_EV = 27.211386018
+BOHR_TO_ANG = 0.529177210903
+
+
+def read_text(path: Path) -> str:
+    return path.read_text(errors="ignore") if path.exists() else ""
+
+
+def detect_backend(path: Path) -> str:
+    root = path if path.is_dir() else path.parent
+    names = {item.name for item in root.iterdir()} if root.is_dir() else set()
+    if "OUTCAR" in names or "POSCAR" in names:
+        return "vasp"
+    if any(root.glob("*.in")) or any(root.glob("*.out")):
+        return "qe"
+    if any(root.glob("*.abi")) or any(root.glob("*.abo")):
+        return "abinit"
+    raise SystemExit(f"Could not detect backend from {path}")
+
+
+def read_energy(path: Path) -> tuple[str, float]:
+    backend = detect_backend(path)
+    root = path if path.is_dir() else path.parent
+    if backend == "vasp":
+        text = read_text(root / "OUTCAR")
+        matches = re.findall(r"TOTEN\s*=\s*([\-0-9.Ee+]+)", text)
+        if not matches:
+            raise SystemExit(f"No TOTEN found in {root / 'OUTCAR'}")
+        return backend, float(matches[-1])
+    if backend == "qe":
+        out_files = sorted(root.glob("*.out"))
+        if not out_files:
+            raise SystemExit(f"No QE output file found in {root}")
+        text = read_text(out_files[0])
+        matches = re.findall(r"!\s+total energy\s+=\s+([\-0-9.DdEe+]+)\s+Ry", text)
+        if not matches:
+            raise SystemExit(f"No total energy found in {out_files[0]}")
+        return backend, float(matches[-1].replace("D", "e").replace("d", "e")) * RY_TO_EV
+    abo_files = sorted(root.glob("*.abo"))
+    out_file = abo_files[0] if abo_files else root / "run.abo"
+    text = read_text(out_file)
+    matches = re.findall(r"\betotal\s+([\-0-9.DdEe+]+)", text)
+    if not matches:
+        raise SystemExit(f"No etotal found in {out_file}")
+    return backend, float(matches[-1].replace("D", "e").replace("d", "e")) * HARTREE_TO_EV
+
+
+def volume(cell: list[list[float]]) -> float:
+    a, b, c = cell
+    return (
+        a[0] * (b[1] * c[2] - b[2] * c[1])
+        - a[1] * (b[0] * c[2] - b[2] * c[0])
+        + a[2] * (b[0] * c[1] - b[1] * c[0])
+    )
+
+
+def read_structure(path: Path) -> tuple[str, dict[str, object]]:
+    backend = detect_backend(path)
+    root = path if path.is_dir() else path.parent
+    if backend == "vasp":
+        poscar = root / "POSCAR"
+        lines = read_text(poscar).splitlines()
+        scale = float(lines[1].split()[0])
+        lattice = [[float(x) * scale for x in lines[i].split()] for i in range(2, 5)]
+        species = lines[5].split()
+        counts = [int(x) for x in lines[6].split()]
+        return backend, {"lattice": lattice, "species": species, "counts": counts, "natoms": sum(counts)}
+    if backend == "qe":
+        in_files = sorted(root.glob("*.in"))
+        if not in_files:
+            raise SystemExit(f"No QE input file found in {root}")
+        text = read_text(in_files[0])
+        species = []
+        counts: dict[str, int] = {}
+        lattice: list[list[float]] = []
+        lines = text.splitlines()
+        for idx, line in enumerate(lines):
+            upper = line.strip().upper()
+            if upper.startswith("ATOMIC_SPECIES"):
+                j = idx + 1
+                while j < len(lines) and lines[j].strip():
+                    label = lines[j].split()[0]
+                    species.append(label)
+                    j += 1
+            if upper.startswith("CELL_PARAMETERS"):
+                lattice = [[float(x) for x in lines[idx + j + 1].split()[:3]] for j in range(3)]
+            if upper.startswith("ATOMIC_POSITIONS"):
+                j = idx + 1
+                while j < len(lines) and lines[j].strip():
+                    label = lines[j].split()[0]
+                    counts[label] = counts.get(label, 0) + 1
+                    j += 1
+                break
+        ordered_counts = [counts.get(label, 0) for label in species]
+        return backend, {"lattice": lattice, "species": species, "counts": ordered_counts, "natoms": sum(ordered_counts)}
+    abi_files = sorted(root.glob("*.abi"))
+    if not abi_files:
+        raise SystemExit(f"No ABINIT input file found in {root}")
+    text = read_text(abi_files[0])
+    acell_match = re.search(r"acell\s+3\*([0-9.]+)", text)
+    rprim_match = re.search(r"rprim\s+(.+?)xred", text, re.DOTALL)
+    if not acell_match or not rprim_match:
+        raise SystemExit(f"Could not parse acell/rprim in {abi_files[0]}")
+    scale = float(acell_match.group(1)) * BOHR_TO_ANG
+    raw = [line.split()[:3] for line in rprim_match.group(1).splitlines() if line.split()]
+    lattice = [[float(x) * scale for x in row] for row in raw[:3]]
+    znucl_match = re.search(r"^\s*znucl\s+(.+)$", text, re.MULTILINE)
+    typat_match = re.search(r"^\s*typat\s+(.+)$", text, re.MULTILINE)
+    species = []
+    counts = []
+    if znucl_match and typat_match:
+        z_values = [int(float(x)) for x in znucl_match.group(1).split()]
+        symbol_map = {3: "Li", 8: "O", 26: "Fe"}
+        species = [symbol_map.get(z, str(z)) for z in z_values]
+        typats = [int(x) for x in typat_match.group(1).split()]
+        for idx, _ in enumerate(species, start=1):
+            counts.append(sum(1 for value in typats if value == idx))
+    return backend, {"lattice": lattice, "species": species, "counts": counts, "natoms": sum(counts)}
diff --git a/scripts/run_regression.py b/scripts/run_regression.py

Original file line number	Diff line number	Diff line change
`@@ -0,0 +1,2 @@`
	`1`	`+ etotal -0.77950703`
	`2`	`+ Calculation completed.`
Original file line number	Diff line number	Diff line change
`@@ -0,0 +1,2 @@`
	`1`	`+ etotal -1.000000`
	`2`	`+ Calculation completed.`