Add multi-species defect chemical potentials

Author: quanshengwu

README.md

@@ -2,7 +2,7 @@
 
 [![CI](https://img.shields.io/github/actions/workflow/status/chatmaterials/defect-analysis/ci.yml?branch=main&label=CI)](https://github.com/chatmaterials/defect-analysis/actions/workflows/ci.yml) [![Release](https://img.shields.io/github/v/release/chatmaterials/defect-analysis?display_name=tag)](https://github.com/chatmaterials/defect-analysis/releases)
 
-Standalone skill for defect-focused DFT result analysis, including automatic defect-type inference and multi-candidate screening.
+Standalone skill for defect-focused DFT result analysis, including automatic defect-type inference, multi-species chemical-potential support, and multi-candidate screening.
 
 Supports VASP, QE, and ABINIT-style pristine/defect inputs.
 
@@ -20,8 +20,9 @@ npx skills add . --list
 python3 scripts/analyze_defect_formation.py fixtures/pristine fixtures/defect --mu -4.0 --temperature-k 1000 --site-density-cm3 1e22 --json
 python3 scripts/analyze_defect_formation.py fixtures/qe/pristine fixtures/qe/defect --mu -4.0 --json
 python3 scripts/analyze_defect_formation.py fixtures/abinit/pristine fixtures/abinit/defect --mu -4.0 --json
+python3 scripts/analyze_defect_formation.py fixtures/substitutional/pristine fixtures/substitutional/defect --mu-term Fe=-6.0 --mu-term Li=-1.5 --json
 python3 scripts/analyze_defect_structure.py fixtures/pristine/POSCAR fixtures/defect/POSCAR --json
-python3 scripts/compare_defect_candidates.py fixtures fixtures/candidates/high-energy-vacancy --mu -4.0 --max-volume-change-percent 5.0 --json
+python3 scripts/compare_defect_candidates.py fixtures fixtures/candidates/high-energy-vacancy --mu -4.0 --max-volume-change-percent 5.0 --temperature-k 1000 --site-density-cm3 1e22 --target-defect-type vacancy-like --target-concentration-cm3 1e10 --json
 python3 scripts/export_defect_report.py fixtures/pristine fixtures/defect --mu -4.0 --temperature-k 1000
 python3 scripts/run_regression.py
 ```

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, automatic defect-type inference, stoichiometric changes, structural distortion around a defect, candidate ranking, and compact markdown reports from pristine and defect calculations. Supports VASP, QE, and ABINIT-style inputs."
+description: "Use when the task is to analyze defect-related DFT results, including neutral defect formation energies, automatic defect-type inference, multi-species chemical-potential terms, stoichiometric changes, abundance estimates, candidate ranking, and compact markdown reports from pristine and defect calculations. Supports VASP, QE, and ABINIT-style inputs."
 ---
 
 # Defect Analysis
@@ -11,8 +11,10 @@ Use this skill for defect-focused post-processing rather than generic workflow s
 
 - estimate a neutral defect formation energy from pristine and defect calculations
 - infer whether the defect looks vacancy-, interstitial-, substitutional-, or complex-like
+- apply multi-species chemical-potential terms for substitutional or complex defects
 - summarize how stoichiometry changes between pristine and defect cells
 - quantify structural or volume change caused by a defect
+- estimate compact abundance metrics under a temperature and site-density assumption
 - rank multiple defect candidates with a compact formation-plus-strain heuristic
 - write a compact defect-analysis report from finished calculations
 
@@ -25,11 +27,11 @@ Supported backends:
 ## Use the bundled helpers
 
 - `scripts/analyze_defect_formation.py`
-  Estimate a neutral defect formation energy, infer the defect type, and optionally estimate an equilibrium fraction.
+  Estimate a neutral defect formation energy, infer the defect type, and optionally estimate abundance metrics.
 - `scripts/analyze_defect_structure.py`
   Compare pristine and defect structures and summarize size, stoichiometry, and inferred defect type.
 - `scripts/compare_defect_candidates.py`
-  Rank multiple defect candidates with a compact formation-plus-strain heuristic.
+  Rank multiple defect candidates with a compact formation-plus-strain-and-abundance heuristic.
 - `scripts/export_defect_report.py`
   Export a markdown defect-analysis report.
 

fixtures/substitutional/defect/OUTCAR

@@ -0,0 +1,2 @@
+ free  energy   TOTEN  =      -19.000000 eV
+ General timing and accounting informations for this job

fixtures/substitutional/defect/POSCAR

@@ -0,0 +1,11 @@
+Li on Fe
+1.0
+5.020000 0.000000 0.000000
+0.000000 5.020000 0.000000
+0.000000 0.000000 5.020000
+Li O
+1 2
+Direct
+0.000000 0.000000 0.000000
+0.500000 0.500000 0.000000
+0.500000 0.000000 0.500000

fixtures/substitutional/pristine/OUTCAR

@@ -0,0 +1,2 @@
+ free  energy   TOTEN  =      -25.000000 eV
+ General timing and accounting informations for this job

fixtures/substitutional/pristine/POSCAR

@@ -0,0 +1,11 @@
+Pristine
+1.0
+5.000000 0.000000 0.000000
+0.000000 5.000000 0.000000
+0.000000 0.000000 5.000000
+Fe O
+1 2
+Direct
+0.000000 0.000000 0.000000
+0.500000 0.500000 0.000000
+0.500000 0.000000 0.500000

references/defect.md

@@ -3,3 +3,4 @@
 - Automatic defect-type inference is a compact structural descriptor, not a substitute for a full site-resolved defect analysis.
 - Neutral formation energies are useful for screening, but charged-defect corrections and Fermi-level dependence can dominate the final physics.
 - Candidate ranking by formation energy and strain is a triage step, not a full defect thermodynamics workflow.
+- Multi-species chemical-potential terms are essential for substitutional defects; a single chemical potential is generally only meaningful for one-species changes.

scripts/analyze_defect_formation.py

@@ -9,14 +9,41 @@
 from defect_io import equilibrium_fraction, infer_species_and_delta, read_energy
 
 
+def parse_mu_terms(raw_terms: list[str] | None) -> dict[str, float]:
+    terms: dict[str, float] = {}
+    for raw in raw_terms or []:
+        if "=" not in raw:
+            raise SystemExit(f"Invalid --mu-term value `{raw}`; expected SPECIES=VALUE")
+        label, value = raw.split("=", 1)
+        terms[label.strip()] = float(value)
+    return terms
+
+
+def classify_abundance(fraction: float | None, concentration: float | None) -> str | None:
+    if concentration is not None:
+        if concentration >= 1e18:
+            return "abundant-like"
+        if concentration >= 1e12:
+            return "dilute-but-relevant"
+        return "trace-like"
+    if fraction is not None:
+        if fraction >= 1e-4:
+            return "abundant-like"
+        if fraction >= 1e-10:
+            return "dilute-but-relevant"
+        return "trace-like"
+    return None
+
+
 def analyze(
     pristine: Path,
     defect: Path,
     species: str | None,
     delta: int | None,
-    mu: float,
+    mu: float | None,
     temperature_K: float | None = None,
     site_density_cm3: float | None = None,
+    mu_terms: dict[str, float] | None = None,
 ) -> dict[str, object]:
     backend_pristine, e_bulk = read_energy(pristine)
     backend_defect, e_def = read_energy(defect)
@@ -25,9 +52,18 @@ def analyze(
     species, delta, defect_type, delta_map = infer_species_and_delta(pristine, defect, species, delta)
     if species is None or delta is None:
         raise SystemExit("Could not infer species and delta from the pristine/defect pair; provide --species and --delta explicitly")
-    e_form = e_def - e_bulk - delta * mu
+    required_species = [label for label, value in delta_map.items() if value != 0]
+    chemical_terms = dict(mu_terms or {})
+    if mu is not None and species not in chemical_terms:
+        chemical_terms[species] = mu
+    missing = [label for label in required_species if label not in chemical_terms]
+    if missing:
+        raise SystemExit(f"Missing chemical potentials for species: {', '.join(missing)}")
+    chemical_sum = sum(delta_map[label] * chemical_terms[label] for label in required_species)
+    e_form = e_def - e_bulk - chemical_sum
     fraction = equilibrium_fraction(e_form, temperature_K) if temperature_K is not None else None
     concentration = fraction * site_density_cm3 if fraction is not None and site_density_cm3 is not None else None
+    abundance_class = classify_abundance(fraction, concentration)
     observations = ["Neutral defect formation energy estimated from pristine and defect total energies."]
     observations.append(f"Defect type inferred as `{defect_type}`.")
     return {
@@ -38,14 +74,16 @@ def analyze(
         "delta_species": delta,
         "species_delta_map": delta_map,
         "defect_type": defect_type,
-        "chemical_potential_eV": mu,
+        "chemical_potential_eV": chemical_terms.get(species),
+        "chemical_potential_terms": chemical_terms,
         "pristine_energy_eV": e_bulk,
         "defect_energy_eV": e_def,
         "formation_energy_eV": e_form,
         "temperature_K": temperature_K,
         "site_density_cm3": site_density_cm3,
         "equilibrium_fraction": fraction,
         "equilibrium_concentration_cm3": concentration,
+        "abundance_class": abundance_class,
         "observations": observations,
     }
 
@@ -56,11 +94,14 @@ def main() -> None:
     parser.add_argument("defect")
     parser.add_argument("--species")
     parser.add_argument("--delta", type=int, help="Change in species count: defect minus pristine.")
-    parser.add_argument("--mu", type=float, required=True, help="Chemical potential for the species in eV.")
+    parser.add_argument("--mu", type=float, help="Chemical potential for the primary species in eV.")
+    parser.add_argument("--mu-term", action="append", help="Chemical potential term in the form SPECIES=VALUE. Repeat as needed.")
     parser.add_argument("--temperature-k", type=float)
     parser.add_argument("--site-density-cm3", type=float)
     parser.add_argument("--json", action="store_true")
     args = parser.parse_args()
+    if args.mu is None and not args.mu_term:
+        raise SystemExit("Provide either --mu or at least one --mu-term")
     payload = analyze(
         Path(args.pristine).expanduser().resolve(),
         Path(args.defect).expanduser().resolve(),
@@ -69,6 +110,7 @@ def main() -> None:
         args.mu,
         args.temperature_k,
         args.site_density_cm3,
+        parse_mu_terms(args.mu_term),
     )
     if args.json:
         print(json.dumps(payload, indent=2))

scripts/compare_defect_candidates.py

@@ -6,51 +6,74 @@
 import json
 from pathlib import Path
 
-from analyze_defect_formation import analyze as analyze_formation
+from analyze_defect_formation import analyze as analyze_formation, parse_mu_terms
 from analyze_defect_structure import analyze as analyze_structure
 
 
 def analyze_case(
     root: Path,
-    mu: float,
+    mu: float | None,
     species: str | None,
     delta: int | None,
     max_volume_change_percent: float,
+    mu_terms: dict[str, float] | None,
+    temperature_K: float | None,
+    site_density_cm3: float | None,
+    target_defect_type: str | None,
+    target_concentration_cm3: float | None,
 ) -> dict[str, object]:
     pristine = root / "pristine"
     defect = root / "defect"
-    formation = analyze_formation(pristine, defect, species, delta, mu)
+    formation = analyze_formation(pristine, defect, species, delta, mu, temperature_K, site_density_cm3, mu_terms)
     structure = analyze_structure(pristine, defect)
     formation_penalty = max(0.0, float(formation["formation_energy_eV"]))
     strain_penalty = max(0.0, abs(float(structure["relative_volume_change_percent"])) - max_volume_change_percent) / 5.0
-    score = formation_penalty + strain_penalty
+    type_penalty = 0.0 if target_defect_type is None or formation["defect_type"].startswith(target_defect_type) else 2.0
+    concentration = formation["equilibrium_concentration_cm3"]
+    abundance_penalty = 0.0
+    if target_concentration_cm3 is not None:
+        abundance_penalty = max(0.0, target_concentration_cm3 - (float(concentration) if concentration is not None else 0.0)) / target_concentration_cm3
+    score = formation_penalty + strain_penalty + type_penalty + abundance_penalty
     return {
         "case": root.name,
         "path": str(root),
         "defect_type": formation["defect_type"],
         "species": formation["species"],
         "delta_species": formation["delta_species"],
         "formation_energy_eV": formation["formation_energy_eV"],
+        "equilibrium_concentration_cm3": concentration,
+        "abundance_class": formation["abundance_class"],
         "relative_volume_change_percent": structure["relative_volume_change_percent"],
         "formation_penalty": formation_penalty,
         "strain_penalty": strain_penalty,
+        "type_penalty": type_penalty,
+        "abundance_penalty": abundance_penalty,
         "screening_score": score,
     }
 
 
 def analyze_cases(
     roots: list[Path],
-    mu: float,
+    mu: float | None,
     species: str | None,
     delta: int | None,
     max_volume_change_percent: float,
+    mu_terms: dict[str, float] | None,
+    temperature_K: float | None,
+    site_density_cm3: float | None,
+    target_defect_type: str | None,
+    target_concentration_cm3: float | None,
 ) -> dict[str, object]:
-    cases = [analyze_case(root, mu, species, delta, max_volume_change_percent) for root in roots]
+    cases = [
+        analyze_case(root, mu, species, delta, max_volume_change_percent, mu_terms, temperature_K, site_density_cm3, target_defect_type, target_concentration_cm3)
+        for root in roots
+    ]
     ranked = sorted(cases, key=lambda item: item["screening_score"])
     return {
         "chemical_potential_eV": mu,
+        "chemical_potential_terms": mu_terms,
         "max_volume_change_percent": max_volume_change_percent,
-        "ranking_basis": "screening_score = formation_penalty + strain_penalty",
+        "ranking_basis": "screening_score = formation_penalty + strain_penalty + type_penalty + abundance_penalty",
         "cases": ranked,
         "best_case": ranked[0]["case"] if ranked else None,
         "observations": [
@@ -62,18 +85,30 @@ def analyze_cases(
 def main() -> None:
     parser = argparse.ArgumentParser(description="Rank defect candidates with a compact formation-plus-strain heuristic.")
     parser.add_argument("paths", nargs="+")
-    parser.add_argument("--mu", type=float, required=True)
+    parser.add_argument("--mu", type=float)
+    parser.add_argument("--mu-term", action="append")
     parser.add_argument("--species")
     parser.add_argument("--delta", type=int)
     parser.add_argument("--max-volume-change-percent", type=float, default=5.0)
+    parser.add_argument("--temperature-k", type=float)
+    parser.add_argument("--site-density-cm3", type=float)
+    parser.add_argument("--target-defect-type")
+    parser.add_argument("--target-concentration-cm3", type=float)
     parser.add_argument("--json", action="store_true")
     args = parser.parse_args()
+    if args.mu is None and not args.mu_term:
+        raise SystemExit("Provide either --mu or at least one --mu-term")
     payload = analyze_cases(
         [Path(path).expanduser().resolve() for path in args.paths],
         args.mu,
         args.species,
         args.delta,
         args.max_volume_change_percent,
+        parse_mu_terms(args.mu_term),
+        args.temperature_k,
+        args.site_density_cm3,
+        args.target_defect_type,
+        args.target_concentration_cm3,
     )
     if args.json:
         print(json.dumps(payload, indent=2))

scripts/export_defect_report.py

@@ -12,10 +12,13 @@
 def screening_note(formation: dict[str, object], structure: dict[str, object]) -> str:
     e_form = float(formation["formation_energy_eV"])
     strain = abs(float(structure["relative_volume_change_percent"]))
+    abundance = formation.get("abundance_class")
     if e_form <= 1.0 and strain <= 5.0:
         return "This defect looks relatively accessible in a compact neutral-defect screen: low formation energy and limited structural swelling."
     if e_form > 2.5:
         return "The neutral formation energy is high enough that this defect is unlikely to be abundant without strong thermodynamic driving forces."
+    if abundance == "trace-like":
+        return "This defect is thermodynamically allowed in the compact model, but the estimated abundance remains trace-like under the supplied conditions."
     if strain > 10.0:
         return "The structural response is large enough that local strain accommodation may be important."
     return "The defect is intermediate in this compact screen and may need charged-defect or concentration analysis before stronger claims."
@@ -29,9 +32,11 @@ def render_markdown(formation: dict[str, object], structure: dict[str, object])
         f"- Defect type: `{formation['defect_type']}`",
         f"- Species: `{formation['species']}`",
         f"- Delta species: `{formation['delta_species']}`",
-        f"- Chemical potential (eV): `{formation['chemical_potential_eV']:.4f}`",
+        f"- Chemical potential terms: `{formation['chemical_potential_terms']}`",
         f"- Formation energy (eV): `{formation['formation_energy_eV']:.4f}`",
         f"- Equilibrium fraction: `{formation['equilibrium_fraction']:.4e}`" if formation["equilibrium_fraction"] is not None else "- Equilibrium fraction: `n/a`",
+        f"- Equilibrium concentration (cm^-3): `{formation['equilibrium_concentration_cm3']:.4e}`" if formation["equilibrium_concentration_cm3"] is not None else "- Equilibrium concentration (cm^-3): `n/a`",
+        f"- Abundance class: `{formation['abundance_class']}`" if formation["abundance_class"] is not None else "- Abundance class: `n/a`",
         "",
         "## Structural Change",
         f"- Pristine atoms: `{structure['natoms_pristine']}`",