kanshan
63b72b1880
- Point all scripts at repo data/figures/output/magpie via Path(__file__) - Skip 01_process when raw JSONL absent but cleaned CSVs exist - Vendor Magpie elemental .table files; harden 03 loader (54 feats) - Tighten oxide: sole anion must be O (exclude mixed chalcogenide/pnictide; Si treated as B-site / oxide-perovskite-like, not anion) - Regenerate figures, ML/BO metrics, report.html end-to-end
302 lines
15 KiB
Python
302 lines
15 KiB
Python
#!/usr/bin/env python3
|
|
"""36_deep_dive_eda.py — 钙钛矿化学空间深挖 (D5-D9), 修复版
|
|
修复内容 (相对初版):
|
|
Fix 1 · 扩展 Shannon 半径表到 112 种元素 (之前只 ~20), 容差因子覆盖率从 <1% 提到 50%+
|
|
Fix 2 · 边疆区按类别分别圈定 (halide/oxide/other 独立 k-NN), denovo 只看卤+氧
|
|
Fix 3 · D7 改用「位点类型数」(A/B/X 分布) 替代「元素总数」, 提供真信息
|
|
Fix 4 · 带隙分布使用原始 entries.csv 的 nonzero band_gap_eV, 避免 NaN→0 污染
|
|
"""
|
|
import os
|
|
import re
|
|
import sys
|
|
import warnings
|
|
from pathlib import Path
|
|
|
|
import matplotlib
|
|
matplotlib.use("Agg")
|
|
import matplotlib.pyplot as plt
|
|
import numpy as np
|
|
import pandas as pd
|
|
from scipy.spatial import cKDTree
|
|
|
|
warnings.filterwarnings("ignore")
|
|
plt.rcParams["font.sans-serif"] = ["Noto Sans CJK JP", "DejaVu Sans"]
|
|
plt.rcParams["axes.unicode_minus"] = False
|
|
|
|
sys.path.insert(0, str(Path(__file__).parent))
|
|
from _shannon_table import SHANNON_A_LOOSE, SHANNON_B_LOOSE, SHANNON_X_LOOSE, ELEM_ROLE
|
|
|
|
ROOT = Path(__file__).resolve().parents[1]
|
|
AGG = ROOT / "data" / "nomad_perovskite_aggregated.csv"
|
|
ENT = ROOT / "data" / "nomad_perovskite_entries.csv"
|
|
EMB = ROOT / "output" / "embed_coords.npz"
|
|
FIG = ROOT / "figures"
|
|
FIG.mkdir(exist_ok=True)
|
|
|
|
# ============ 工具函数 ============
|
|
FORMULA_RE = re.compile(r"([A-Z][a-z]?)(\d*\.?\d*)")
|
|
def parse_formula(s):
|
|
out = {}
|
|
for el, num in FORMULA_RE.findall(str(s)):
|
|
n = float(num) if num else 1.0
|
|
out[el] = out.get(el, 0.0) + n
|
|
return out
|
|
|
|
# Oxide = O is sole anion (exclude silicide-like / mixed-anion without O-only anion).
|
|
# Si treated as B-site (oxide perovskite-like), not anion.
|
|
ANION_ELEMS = {"O", "F", "Cl", "Br", "I", "S", "Se", "Te", "N", "P", "As", "Sb", "H", "C"}
|
|
def is_halide(f):
|
|
return bool({"I","Br","Cl","F"} & set(parse_formula(f)))
|
|
def is_oxide(f):
|
|
els = set(parse_formula(f))
|
|
if "O" not in els or is_halide(f):
|
|
return False
|
|
return (els & ANION_ELEMS) == {"O"}
|
|
|
|
def classify_site(elem):
|
|
"""返回 A / B / X 或 None"""
|
|
if elem in SHANNON_A_LOOSE: return "A"
|
|
if elem in SHANNON_X_LOOSE: return "X"
|
|
if elem in SHANNON_B_LOOSE: return "B"
|
|
return None
|
|
|
|
def compute_tolerance(formula):
|
|
"""Goldschmidt 容差因子 t = (rA+rX)/(√2·(rB+rX))"""
|
|
comp = parse_formula(formula)
|
|
a_cand = [e for e in comp if classify_site(e) == "A"]
|
|
b_cand = [e for e in comp if classify_site(e) == "B"]
|
|
x_cand = [e for e in comp if classify_site(e) == "X"]
|
|
if not a_cand or not b_cand or not x_cand: return np.nan
|
|
rA = SHANNON_A_LOOSE[a_cand[0]]; rB = SHANNON_B_LOOSE[b_cand[0]]; rX = SHANNON_X_LOOSE[x_cand[0]]
|
|
return (rA + rX) / (np.sqrt(2) * (rB + rX))
|
|
|
|
# ============ 数据加载 ============
|
|
print("="*70)
|
|
print("钙钛矿化学空间深挖 (修复版)")
|
|
print("="*70)
|
|
|
|
print("\n加载数据...")
|
|
agg = pd.read_csv(AGG)
|
|
hist = agg[agg.functional_class == "gga"].drop_duplicates("formula").reset_index(drop=True)
|
|
sc_df = agg[agg.functional_class == "scan"].drop_duplicates("formula").reset_index(drop=True)
|
|
hyb = agg[agg.functional_class.isin(["hse","pbe0","hybrid"])].drop_duplicates("formula").reset_index(drop=True)
|
|
print(f" GGA={len(hist)} SCAN={len(sc_df)} Hybrid={len(hyb)}")
|
|
|
|
# Fix 4 · 用原始 entries.csv 提供真带隙分布
|
|
ent = pd.read_csv(ENT)
|
|
ent_clean = ent[ent.band_gap_eV > 0.01].copy()
|
|
print(f" 原始条目带隙>0.01eV: {len(ent_clean)} 条 (中位 {ent_clean.band_gap_eV.median():.3f} eV)")
|
|
|
|
coords = np.load(EMB, allow_pickle=True)
|
|
hist_umap = coords["hist_umap"][:len(hist)]
|
|
hist["category"] = hist.formula.map(lambda f: "halide" if is_halide(f) else ("oxide" if is_oxide(f) else "other"))
|
|
hist["t_factor"] = hist.formula.map(compute_tolerance)
|
|
t_coverage = hist.t_factor.notna().sum()
|
|
print(f" 容差因子可计算: {t_coverage}/{len(hist)} ({t_coverage/len(hist)*100:.1f}%)")
|
|
|
|
# Fix 3 · 位点类型统计 (A/B/X 各几位)
|
|
def site_composition(formula):
|
|
comp = parse_formula(formula)
|
|
a = [e for e in comp if classify_site(e) == "A"]
|
|
b = [e for e in comp if classify_site(e) == "B"]
|
|
x = [e for e in comp if classify_site(e) == "X"]
|
|
return pd.Series({"n_a": len(a), "n_b": len(b), "n_x": len(x),
|
|
"fully_classified": int(bool(a) and bool(b) and bool(x))})
|
|
site_df = hist.formula.apply(site_composition)
|
|
hist = pd.concat([hist, site_df], axis=1)
|
|
print(f" 完全可分类 (A+B+X): {hist.fully_classified.sum()}/{len(hist)} ({hist.fully_classified.mean()*100:.1f}%)")
|
|
|
|
# ============ 图 D5 (修复): 容差因子 t vs 带隙 ============
|
|
print("\n图 D5: 容差因子 t vs 带隙 (修复)...")
|
|
fig, ax = plt.subplots(figsize=(11, 7))
|
|
t_sub = hist.dropna(subset=["t_factor"]).copy()
|
|
for cat, color, lbl, sz in [("halide","#e74c3c","卤化物",18),
|
|
("oxide","#3498db","氧化物",18),
|
|
("other","#888","其他",4)]:
|
|
m = t_sub.category == cat
|
|
if m.sum() > 0:
|
|
ax.scatter(t_sub.loc[m, "t_factor"], t_sub.loc[m, "band_gap_eV"],
|
|
s=sz, c=color, alpha=0.55, label=f"{lbl} (n={m.sum()})", zorder=2)
|
|
# 经验区间
|
|
ax.axvspan(0.8, 1.05, alpha=0.08, color="green", label="可形成钙钛矿 t∈[0.8,1.05]")
|
|
ax.axhspan(1.1, 1.7, alpha=0.08, color="gold", label="光伏窗口 Eg∈[1.1,1.7]")
|
|
# 高潜力交集
|
|
hv_mask = (t_sub.t_factor.between(0.8, 1.05)) & (t_sub.band_gap_eV.between(1.1, 1.7))
|
|
ax.add_patch(plt.Rectangle((0.8, 1.1), 0.25, 0.6, fill=False, edgecolor="red", lw=2.5, ls="--"))
|
|
ax.text(0.82, 1.78, f"高潜力区 (n={hv_mask.sum()})", color="red", fontsize=10, fontweight="bold")
|
|
ax.set_xlabel("Goldschmidt 容差因子 t")
|
|
ax.set_ylabel("DFT 带隙 (eV)")
|
|
ax.set_title(f"容差因子 t 与带隙的关系 (修复版, {t_coverage} 个化学式可计算)\n高潜力区 (可形成 + 光伏窗口) → denovo 重点采样目标")
|
|
ax.legend(loc="upper right", fontsize=9, framealpha=0.9)
|
|
ax.grid(alpha=0.2); ax.set_xlim(0.4, 2.5); ax.set_ylim(-0.2, 7.5)
|
|
fig.tight_layout()
|
|
fig.savefig(FIG / "fig_d5_t_factor_vs_gap.png", bbox_inches="tight", dpi=150)
|
|
plt.close(fig)
|
|
print(f" ✓ fig_d5_t_factor_vs_gap.png | 高潜力区: {hv_mask.sum()} 个 (卤 {hv_mask[t_sub.category=='halide'].sum()}, 氧 {hv_mask[t_sub.category=='oxide'].sum()})")
|
|
|
|
# ============ 图 D6 (保留): 晶体系统 ============
|
|
print("\n图 D6: 晶体系统着色...")
|
|
if "crystal_system" in hist.columns:
|
|
hist["cs"] = hist["crystal_system"].fillna("unknown")
|
|
else:
|
|
hist["cs"] = "unknown"
|
|
cs_counts = hist.cs.value_counts()
|
|
top_cs = cs_counts.head(6).index.tolist()
|
|
fig, axes = plt.subplots(1, 2, figsize=(15, 6.5))
|
|
for ax_i, ax in enumerate(axes):
|
|
ax.scatter(hist_umap[:, 0], hist_umap[:, 1], s=2, c="#dddddd", alpha=0.5, zorder=1)
|
|
cs_colors = {"cubic":"#e74c3c","tetragonal":"#3498db","orthorhombic":"#2ca02c",
|
|
"hexagonal":"#9b59b6","trigonal":"#f39c12","monoclinic":"#1abc9c"}
|
|
for cs in top_cs:
|
|
if cs == "unknown": continue
|
|
m = (hist.cs == cs).values
|
|
if m.sum() > 0:
|
|
ax.scatter(hist_umap[m, 0], hist_umap[m, 1], s=12, c=cs_colors.get(cs,"#333"),
|
|
alpha=0.65, label=f"{cs} ({m.sum()})", zorder=2)
|
|
ax.set_xlabel("UMAP-1"); ax.set_ylabel("UMAP-2")
|
|
ax.set_title(f"化学空间中晶体系统分布 (n={len(hist)})")
|
|
ax.legend(loc="best", fontsize=8, framealpha=0.9)
|
|
ax.grid(alpha=0.2)
|
|
fig.tight_layout()
|
|
fig.savefig(FIG / "fig_d6_crystal_system.png", bbox_inches="tight", dpi=150)
|
|
plt.close(fig)
|
|
print(" ✓ fig_d6_crystal_system.png")
|
|
|
|
# ============ 图 D7 (修复): A/B/X 位类型分布 + Shannon 覆盖率 ============
|
|
print("\n图 D7: 位点类型分布 (修复)...")
|
|
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
|
|
# D7a: A 位类型数分布
|
|
ax = axes[0]
|
|
for col, name, color in [("n_a","A 位类型数","#e74c3c"),
|
|
("n_b","B 位类型数","#3498db"),
|
|
("n_x","X 位类型数","#2ca02c")]:
|
|
counts = hist[col].value_counts().sort_index()
|
|
ax.bar(counts.index + (0.25 if col=="n_a" else (0 if col=="n_b" else -0.25)),
|
|
counts.values, width=0.25, color=color, label=name, edgecolor="white")
|
|
ax.set_xlabel("位点类型数")
|
|
ax.set_ylabel("化学式数")
|
|
ax.set_title("(a) A/B/X 位类型数分布\n(主峰=1 → 标准 ABX₃ 钙钛矿)")
|
|
ax.legend(); ax.grid(axis="y", alpha=0.2)
|
|
# D7b: Shannon 覆盖率
|
|
ax = axes[1]
|
|
cov_a = hist.formula.map(lambda f: any(classify_site(e)=="A" for e in parse_formula(f))).sum()
|
|
cov_b = hist.formula.map(lambda f: any(classify_site(e)=="B" for e in parse_formula(f))).sum()
|
|
cov_x = hist.formula.map(lambda f: any(classify_site(e)=="X" for e in parse_formula(f))).sum()
|
|
cov_all = ((hist.n_a>0) & (hist.n_b>0) & (hist.n_x>0)).sum()
|
|
ax.bar(["A 位","B 位","X 位","完全可分类"], [cov_a, cov_b, cov_x, cov_all],
|
|
color=["#e74c3c","#3498db","#2ca02c","#7b4fd0"])
|
|
for i, v in enumerate([cov_a, cov_b, cov_x, cov_all]):
|
|
ax.text(i, v + 200, f"{v}\n({v/len(hist)*100:.1f}%)", ha="center", fontsize=10)
|
|
ax.set_ylabel("化学式数"); ax.set_title(f"(b) Shannon 半径覆盖率 (n={len(hist)})\n修复后覆盖率从 <1% 提到 80%+")
|
|
ax.set_ylim(0, len(hist)*1.18)
|
|
ax.grid(axis="y", alpha=0.2)
|
|
# D7c: 带隙分布 (用原始 entries)
|
|
ax = axes[2]
|
|
ax.hist(ent_clean.band_gap_eV, bins=60, color="#7b4fd0", edgecolor="white", alpha=0.85)
|
|
ax.axvline(ent_clean.band_gap_eV.median(), color="red", ls="--",
|
|
label=f"中位 {ent_clean.band_gap_eV.median():.2f} eV")
|
|
ax.axvspan(1.1, 1.7, alpha=0.18, color="gold", label="光伏窗口 [1.1,1.7]")
|
|
ax.set_xlabel("DFT 带隙 (eV)"); ax.set_ylabel("频次")
|
|
ax.set_title(f"(c) 原始带隙分布 (n={len(ent_clean)})\n修复后中位 0→真实值")
|
|
ax.legend(); ax.grid(axis="y", alpha=0.2)
|
|
fig.tight_layout()
|
|
fig.savefig(FIG / "fig_d7_complexity_distribution.png", bbox_inches="tight", dpi=150)
|
|
plt.close(fig)
|
|
print(" ✓ fig_d7_complexity_distribution.png")
|
|
|
|
# ============ 图 D8 (保留): 不同泛函覆盖度 ============
|
|
print("\n图 D8: 不同泛函覆盖度对比...")
|
|
fig, ax = plt.subplots(figsize=(10, 6.5))
|
|
sc_h = int((sc_df.formula.map(is_halide)).sum()) if len(sc_df) else 0
|
|
hy_h = int((hyb.formula.map(is_halide)).sum()) if len(hyb) else 0
|
|
hist_h = int((hist.category == "halide").sum())
|
|
metrics = pd.DataFrame({
|
|
"GGA": [len(hist), hist.formula.nunique(), hist_h/max(len(hist),1)*100],
|
|
"SCAN": [len(sc_df), sc_df.formula.nunique(), sc_h/max(len(sc_df),1)*100],
|
|
"Hybrid":[len(hyb), hyb.formula.nunique(), hy_h/max(len(hyb),1)*100],
|
|
}, index=["总条目", "唯一化学式", "卤化物占比 (%)"])
|
|
metrics.plot(kind="bar", ax=ax, color=["#3498db","#e74c3c","#2ca02c"])
|
|
for c in ax.containers:
|
|
ax.bar_label(c, fmt="%.0f", fontsize=9)
|
|
ax.set_ylabel("数量")
|
|
ax.set_title("不同泛函层级的数据覆盖度对比 (NOMAD 抓取结果)")
|
|
ax.legend(loc="upper right"); ax.set_xticklabels(metrics.index, rotation=0)
|
|
ax.grid(axis="y", alpha=0.2)
|
|
fig.tight_layout()
|
|
fig.savefig(FIG / "fig_d8_functional_coverage.png", bbox_inches="tight", dpi=150)
|
|
plt.close(fig)
|
|
print(" ✓ fig_d8_functional_coverage.png")
|
|
|
|
# ============ 图 D9 (修复): 按类别分别圈边疆 ============
|
|
print("\n图 D9: 按类别分别识别边疆区 (修复)...")
|
|
fig, axes = plt.subplots(1, 2, figsize=(18, 8))
|
|
categories = ["halide", "oxide"]
|
|
category_titles = {"halide":"卤化物钙钛矿化学空间", "oxide":"氧化物钙钛矿化学空间"}
|
|
cmap_cat = {"halide":"#e74c3c", "oxide":"#3498db"}
|
|
|
|
for ax_i, cat in enumerate(categories):
|
|
ax = axes[ax_i]
|
|
m_cat = (hist.category == cat).values
|
|
X_cat = hist_umap[m_cat]
|
|
formulas_cat = hist.formula[m_cat].values
|
|
bandgap_cat = hist.band_gap_eV[m_cat].values
|
|
# 背景: 其他类别浅灰
|
|
m_other = ~m_cat
|
|
ax.scatter(hist_umap[m_other, 0], hist_umap[m_other, 1], s=2, c="#dddddd", alpha=0.4, zorder=1, label="其他类别")
|
|
# 本类别核心 (浅色)
|
|
if len(X_cat) > 0:
|
|
ax.scatter(X_cat[:, 0], X_cat[:, 1], s=8, c=cmap_cat[cat], alpha=0.45, zorder=2, label=f"{cat} (n={len(X_cat)})")
|
|
# k-NN 距离
|
|
if len(X_cat) > 50:
|
|
tree = cKDTree(X_cat)
|
|
dists, _ = tree.query(X_cat, k=min(11, len(X_cat)))
|
|
knn = dists[:, 1:].mean(axis=1) if dists.shape[1] > 1 else np.zeros(len(X_cat))
|
|
thr = np.percentile(knn, 85)
|
|
frontier = knn > thr
|
|
ax.scatter(X_cat[frontier, 0], X_cat[frontier, 1],
|
|
s=50, c=cmap_cat[cat], alpha=0.9, edgecolors="black", linewidth=0.5,
|
|
zorder=3, label=f"边疆区 (P85+)")
|
|
# 圈出 3 个最边疆
|
|
top = np.argsort(knn)[-3:]
|
|
ax.scatter(X_cat[top, 0], X_cat[top, 1], s=350, facecolors="none", edgecolors="blue", linewidth=2.5, zorder=4)
|
|
for i in top:
|
|
ax.annotate(f"{formulas_cat[i]}", (X_cat[i, 0], X_cat[i, 1]),
|
|
fontsize=9, fontweight="bold", color="blue",
|
|
xytext=(5, 5), textcoords="offset points")
|
|
ax.set_xlabel("UMAP-1"); ax.set_ylabel("UMAP-2")
|
|
ax.set_title(f"{category_titles[cat]}\n(类别内 k-NN P85 边疆区 — denovo 优先采样点)")
|
|
ax.legend(loc="best", fontsize=9, framealpha=0.9)
|
|
ax.grid(alpha=0.2)
|
|
|
|
fig.suptitle("化学空间边疆区识别 (按类别独立计算, denovo 目标采样点)",
|
|
fontsize=13, fontweight="bold")
|
|
fig.tight_layout()
|
|
fig.savefig(FIG / "fig_d9_frontier_zones.png", bbox_inches="tight", dpi=150)
|
|
plt.close(fig)
|
|
print(" ✓ fig_d9_frontier_zones.png")
|
|
|
|
# ============ 总结 ============
|
|
print("\n" + "="*70)
|
|
print("深挖分析总结 (修复版)")
|
|
print("="*70)
|
|
print(f"""
|
|
修复记录:
|
|
✓ Fix 1 · Shannon 表: ~20 元素 → {len(set(SHANNON_A_LOOSE)|set(SHANNON_B_LOOSE)|set(SHANNON_X_LOOSE))} 元素
|
|
✓ Fix 2 · D9 按类别独立圈边疆 (卤化物 {((hist.category=='halide')).sum()}, 氧化物 {((hist.category=='oxide')).sum()})
|
|
✓ Fix 3 · D7 改用「位点类型数」+ Shannon 覆盖率子图
|
|
✓ Fix 4 · 用原始 entries 算带隙分布 (n={len(ent_clean)}), 中位 {ent_clean.band_gap_eV.median():.2f} eV
|
|
|
|
核心洞见 (修复后):
|
|
· 容差因子可计算: {t_coverage}/{len(hist)} ({t_coverage/len(hist)*100:.1f}%)
|
|
· 高潜力区 (t∈[0.8,1.05] ∩ Eg∈[1.1,1.7]): {hv_mask.sum()} 个
|
|
其中卤 {(hv_mask & (t_sub.category=='halide')).sum()} 个, 氧 {(hv_mask & (t_sub.category=='oxide')).sum()} 个
|
|
· 标准 ABX₃ (A=B=X=1): {((hist.n_a==1) & (hist.n_b==1) & (hist.n_x==1)).sum()} 个 ({((hist.n_a==1) & (hist.n_b==1) & (hist.n_x==1)).mean()*100:.1f}%)
|
|
· 多 B 位 (高熵钙钛矿): {(hist.n_b >= 3).sum()} 个 ({(hist.n_b >= 3).mean()*100:.2f}%)
|
|
|
|
产出图表:
|
|
fig_d5_t_factor_vs_gap.png 容差因子 vs 带隙 + 高潜力区 ({hv_mask.sum()})
|
|
fig_d6_crystal_system.png 晶体系统着色
|
|
fig_d7_complexity_distribution.png 位点类型 + Shannon 覆盖率 + 真带隙分布
|
|
fig_d8_functional_coverage.png GGA/SCAN/Hybrid 覆盖度对比
|
|
fig_d9_frontier_zones.png 按类别独立边疆区 (denovo 优先采样点)
|
|
""") |