pymatgen 实战教程：材料建模、DFT 前后处理与高通量计算

环境：Windows + Miniconda，Python 3.14.4，VS Code，纯 .py 脚本
基准版本：pymatgen 2025.x（当前最新稳定版）
本文所有 API 均基于 2025.x 真实用法编写；凡涉及联网/API key、或不同版本间易变的接口，都会明确标注。建议写作时对照官方文档 pymatgen.org 核对版本差异。

目录

  1. 简介与定位
  2. 安装与版本检查
  3. 核心对象：Element / Composition / Lattice / Structure / Molecule
  4. 读写结构：CIF、POSCAR、多格式互转
  5. 构建与修改结构
  6. 表面与 Slab：催化吸附建模
  7. Materials Project 数据获取：MPRester
  8. PhaseDiagram：相图与凸包稳定性
  9. 与 VASP 对接：Vasprun / Poscar 后处理
  10. 分析工具：对称性、近邻、键长
  11. 科研实战示例
  12. 常见坑
  13. 延伸资源

1. 简介与定位

pymatgen（Python Materials Genomics）是材料科学计算领域事实上的标准 Python 库，由 Materials Project 团队维护。在 DFT 工作流中，它通常扮演三类角色：

  1. 结构建模与操纵：构建、读写、变换晶体结构（POSCAR/CIF），切表面、建超胞、掺杂、构建吸附构型。
  2. DFT 前后处理：解析 VASP 输出（vasprun.xml、OUTCAR、OSZICAR），提取能量、能带、态密度，生成 VASP 输入。
  3. 高通量与数据库：对接 Materials Project 数据库，用 PhaseDiagram 做相图/凸包稳定性分析，批量筛选候选材料。

在科研场景里，pymatgen 最常用的部分是：

  • 用 Structure / SlabGenerator 搭催化剂表面与吸附构型；
  • 用 Vasprun 读能量算吸附能/形成能；
  • 用 PhaseDiagram 筛稳定相；
  • 用 MPRester 从 Materials Project 拉结构、批量拿形成能数据。

pymatgen 底层依赖 numpy / scipy / spglib（对称性）/ matplotlib（绘图，可选）。核心对象大多在 pymatgen.core，IO 在 pymatgen.io.*，分析在 pymatgen.analysis.*。

2. 安装与版本检查

2.1 安装

推荐用 conda（conda-forge 会一并解决 numpy/scipy/spglib 等编译依赖）：

# 新建环境（若 Python 3.14 下 spglib/numpy 尚无预编译 wheel，可退回 3.12/3.13）
conda create -n mat python=3.12 -y
conda activate mat
conda install -c conda-forge pymatgen -y

或用 pip：

pip install pymatgen

提示：Python 3.14 较新，某些二进制依赖（spglib、scipy）可能暂时没有对应 wheel。若 pip install pymatgen 报编译错误，优先用 conda-forge，或改用 3.12/3.13 环境。这属于生态兼容问题，与 pymatgen 本身无关。

2.2 版本检查

# version_check.py
from importlib.metadata import version

print("pymatgen version:", version("pymatgen"))

# pymatgen.core 内也暴露版本号
import pymatgen.core
print("core.__version__:", pymatgen.core.__version__)

建议脚本里对关键版本做一次断言，避免团队环境不一致导致 API 差异：

from importlib.metadata import version
assert version("pymatgen") >= "2024.0", "pymatgen 版本过旧"

2.3 可选：配置 POTCAR 库（仅当你需要让 pymatgen 自动写 POTCAR）

pymatgen 写 VASP 输入集、生成 POTCAR 时，需要知道赝势库路径。可用 PMG_VASP_PSP_DIR 环境变量，或写 ~/.pmgrc.yaml。若只做后处理（读 vasprun.xml），无需配置。

# 以环境变量为例
export PMG_VASP_PSP_DIR="D:/vasp_pot/potpaw_PBE.54"

3. 核心对象

pymatgen 的核心对象都从 pymatgen.core 导入。最常用：Element、Species、Composition、Lattice、Structure、Molecule。

3.1 Element（元素）

from pymatgen.core import Element

fe = Element("Fe")
print(fe.Z)              # 原子序数 26
print(fe.symbol)         # 'Fe'
print(fe.X)              # 电负性（Pauling）
print(fe.atomic_mass)    # 相对原子质量
print(fe.group, fe.row)  # 族 / 周期
print(fe.is_transition_metal)   # 过渡金属判断（催化里常用）

Element 也可带自旋极化相关属性（如 Element("Fe").spin 缺省未定）。单原子催化筛选中，常用 is_transition_metal / 电负性做初筛。

3.2 Species（带氧化态的元素）

当需要区分价态（如 Fe2+ vs Fe3+）时用 Species：

from pymatgen.core import Species

fe3 = Species("Fe", 3)     # Fe3+
fe2 = Species("Fe", 2)     # Fe2+
print(fe3.oxi_state)       # 3
print(str(fe3))            # 'Fe3+'

结构里通常用 Element 即可；需要价态相关分析（如键价和、氧化态标记）时才用 Species。

3.3 Composition（化学式/组分）

from pymatgen.core import Composition

c = Composition("LiFePO4")
print(c.reduced_formula)        # 'LiFe(PO4)1' 或类似的最简式
print(c.weight)                 # 分子量（相对）
print(c.num_atoms)              # 总原子数 = 7

# 从字典构造
c2 = Composition({"Fe": 2, "O": 3})
print(c2.formula)               # 'Fe2 O3'
print(c2.get_atomic_fraction("Fe"))   # 0.4
print(c2.get_wt_fraction("Fe"))       # 质量分数

# 元素计数
print(c.get_el_amt_dict())      # {Element('Li'): 1.0, ...}

Composition 是相图、形成能、MPRester 搜索里的核心类型，也支持直接加减、比较。

3.4 Lattice（晶格）

from pymatgen.core import Lattice

# 立方晶格
l = Lattice.cubic(3.615)   # 面心 Cu 的晶格常数
print(l.volume)            # 体积（Å^3）
print(l.a, l.b, l.c)       # 晶格常数
print(l.alpha, l.beta, l.gamma)   # 角度（度）

# 通用三斜参数
l2 = Lattice.from_parameters(4.0, 4.0, 4.0, 90, 90, 90)

# 直接从 3x3 矩阵（行向量为晶格基矢）
l3 = Lattice([[3.5, 0, 0], [0, 3.5, 0], [0, 0, 3.5]])
print(l3.matrix)           # 基矢矩阵

3.5 Structure（晶体结构）—— 重中之重

Structure 是最核心的对象。构造方式很多，最常用三种：

from pymatgen.core import Structure, Lattice

# 方式 1：晶格 + 元素列表 + 坐标列表（默认分数坐标）
lattice = Lattice.cubic(3.615)
s = Structure(lattice, ["Cu", "Cu"], [[0, 0, 0], [0.5, 0.5, 0.5]])

# 方式 2：每个原子写成 [元素, 坐标] 的联合列表
s2 = Structure(lattice, ["Cu", [0, 0, 0]], ["Cu", [0.5, 0.5, 0.5]])

# 方式 3：用空间群 + 不等价位置自动展开（推荐用于高对称结构）
s3 = Structure.from_spacegroup(
    "Fm-3m",             # 空间群符号或编号
    Lattice.cubic(3.615),
    ["Cu"],
    [[0, 0, 0]],         # 不等价 Wyckoff 位置的分数坐标
)
print(s3.num_sites)      # 4 —— Fm-3m 自动生成 4 个 Cu

访问结构信息：

print(s.num_sites)        # 原子数
print(s.formula)          # 化学式
print(s.composition)      # Composition 对象
print(s.lattice)          # Lattice 对象
print(s.volume)           # 体积
print(s.density)          # 密度 g/cm^3
print(s.cart_coords)      # 笛卡尔坐标数组 (N,3)
print(s.frac_coords)      # 分数坐标数组 (N,3)

# 遍历原子（site）
for site in s:
    print(site.index, site.specie, site.frac_coords, site.coords)

# 按索引取 site
site0 = s[0]
print(site0.specie)              # Element('Cu')
print(site0.coords)              # 笛卡尔坐标
print(site0.frac_coords)         # 分数坐标

不可变版本 IStructure：如果只读、不想被意外修改，可用 pymatgen.core.IStructure（构造方式相同，但不支持增删改）。用于缓存/字典 key 时更安全。

3.6 Molecule（分子）

from pymatgen.core import Molecule

co2 = Molecule(["C", "O", "O"],
               [[0, 0, 0], [-1.16, 0, 0], [1.16, 0, 0]])
print(co2.formula)
print(co2.composition)

co2.center()               # 平移到质心为原点
print(co2.get_centroid())

Molecule 用于吸附质（CO、CO2、H2O 等）建模、分子筛/MOF 中的客体分子等。与 Structure 共享很多接口。

4. 读写结构

pymatgen 的 IO 分散在 pymatgen.io.* 子包：CIF 在 pymatgen.io.cif，VASP 在 pymatgen.io.vasp。

4.1 通用读取器 Structure.from_file

Structure.from_file 会根据扩展名自动识别格式，是最省心的入口：

from pymatgen.core import Structure

# 读 CIF
s1 = Structure.from_file("input.cif")
# 读 POSCAR
s2 = Structure.from_file("POSCAR")
# 读 CONTCAR
s3 = Structure.from_file("CONTCAR")

它内部对 CIF 走 CifParser，对 POSCAR/CONTCAR 走 Poscar。对普通结构文件足够用。

4.2 CIF：CifParser / CifWriter

from pymatgen.io.cif import CifParser, CifWriter

# 解析 CIF（一个 CIF 文件可能含多个结构）
parser = CifParser("example.cif")
structures = parser.parse_structures(primitive=False)   # 返回 list[Structure]
s = structures[0]

# primitive=True 会尝试约化到原胞
prim = parser.parse_structures(primitive=True)[0]

# 写 CIF
CifWriter(s).write_file("out.cif")

# 或直接用 Structure 的方法（按扩展名推断格式）
s.to(filename="out.cif")

4.3 POSCAR：Poscar

from pymatgen.io.vasp import Poscar

s = Structure.from_file("CONTCAR")   # 已有结构

# 写 POSCAR
poscar = Poscar(s)
poscar.write_file("POSCAR")

# 读 POSCAR
p = Poscar.from_file("POSCAR")
s_read = p.structure

Poscar 还支持写入时给每个元素指定速度（分子动力学）、写 selective dynamics 标志等高级选项（Poscar(s, sort_structure=False, site_symbols=...)）。

4.4 多格式互转

from pymatgen.core import Structure
from pymatgen.io.cif import CifWriter

s = Structure.from_file("POSCAR")

# 统一写法：to() 按扩展名决定格式，返回 None（写入文件）
s.to(filename="out.cif")
s.to(filename="out.cssr")

# 若要得到格式字符串（不落盘），用 fmt 参数
poscar_str = s.to(fmt="poscar")   # 返回 POSCAR 文本
cif_str = s.to(fmt="cif")         # 返回 CIF 文本
print(poscar_str[:80])

批量互转的实用脚本：

# convert.py —— 把一个文件夹里的 POSCAR 全转成 CIF
from pathlib import Path
from pymatgen.core import Structure

for f in Path("poscars").glob("POSCAR*"):
    s = Structure.from_file(f)
    s.to(filename=f"cifs/{f.name}.cif")

5. 构建与修改结构

构建催化模型、做掺杂、切超胞时，下面这些方法每天都会用到。

5.1 增删改原子

from pymatgen.core import Structure, Lattice

s = Structure.from_spacegroup("Fm-3m", Lattice.cubic(3.615), ["Cu"], [[0, 0, 0]])

# 添加原子（默认分数坐标）
s.append("O", [0.5, 0.5, 0.5], coords_are_cartesian=False)

# 在指定索引插入
s.insert(0, "O", [0.25, 0.25, 0.25], coords_are_cartesian=False)

# 替换某个位置的元素/坐标
s.replace(0, "N", [0.1, 0.1, 0.1], coords_are_cartesian=False)

# 删除指定索引的原子
s.remove_sites([0, 2])          # 删除第 0、2 号原子
# 删除指定元素的所有原子
s.remove_species(["O"])

print(s.formula)

5.2 平移 / 缩放

# 平移若干原子（分数坐标）
s.translate_sites([0, 1], [0.5, 0, 0], frac_coords=True, to_unit_cell=True)

# 整体缩放晶格体积到指定值
s.scale_lattice(s.volume * 2)    # 体积翻倍

5.3 掺杂 / 替换（substitute vs replace_species）

# 推荐：replace_species，接受 dict 映射，稳定且直观
s_nico = s.replace_species({"Cu": "Ni"})
print(s_nico.formula)

# substitute 同样可用（接受 dict 映射）。不同版本曾调整签名，不确定时以 replace_species 为准
s_nico2 = s.substitute({"Cu": "Co"})

掺杂质子的常见套路：先 make_supercell，再 replace_species 或 replace 某几个 site。

5.4 超胞

# 2x2x2 超胞，两种等价写法
super1 = s.make_supercell([[2, 0, 0], [0, 2, 0], [0, 0, 2]])
super2 = s * (2, 2, 2)

print(super1.num_sites)

5.5 排序 / 扰动 / 排序副本

s.sort()                    # 就地按电负性等排序
s_sorted = s.get_sorted_structure()   # 返回排序后的副本（不修改原对象）

s.perturb(0.01)             # 随机扰动原子坐标 0.01 Å（建初始构型/打破对称常用）

5.6 氧化态

s.add_oxidation_state_by_guess()   # 按电负性猜测氧化态
print([site.specie.oxi_state for site in s])
s.remove_oxidation_states()

6. 表面与 Slab

切表面是催化建模的第一步。核心类是 pymatgen.core.surface.SlabGenerator。

from pymatgen.core import Structure, Lattice
from pymatgen.core.surface import SlabGenerator

# 1. 先建体相
cu = Structure.from_spacegroup("Fm-3m", Lattice.cubic(3.615), ["Cu"], [[0, 0, 0]])

# 2. 切 (111) 表面
slabgen = SlabGenerator(
    cu,
    miller_index=(1, 1, 1),   # 米勒指数
    min_slab_size=10.0,       # 最小 slab 厚度（Å）
    min_vacuum_size=15.0,     # 最小真空层厚度（Å）
    lll_reduce=True,          # 用 LLL 约化晶格，让超胞尽量"方正"
    center_slab=True,         # 把 slab 居中、真空均分在上下
)
slab = slabgen.get_slab(shift=0)      # 取某个 shift 对应的终结面
slabs = slabgen.get_slabs()           # 所有对称不等价的终结面
print("终结面数量:", len(slabs))
print(slab.num_sites, slab.lattice.c)

关键点：

  • min_vacuum_size 至少 12–15 Å，做带电体系/偶极校正时要更大。
  • 表面建模往往还要 make_supercell 扩成 p(x n) 超胞，避免吸附质自相互作用。
  • get_slabs() 返回多个终结面，通常用 get_slab(shift=i) 逐个选取，比较表面能后挑最稳定的。

加吸附质：AdsorbateSiteFinder（高级）

from pymatgen.core import Molecule
from pymatgen.analysis.adsorption import AdsorbateSiteFinder

asf = AdsorbateSiteFinder(slab)
sites = asf.find_adsorption_sites()    # dict，含 ontop/bridge/hollow 等吸附位坐标

co = Molecule(["C", "O"], [[0, 0, 0], [0, 0, 1.15]])   # CO 吸附质
ads_slab = asf.add_adsorbate(co, sites["ontop"][0])    # 放到顶位
ads_slab.to(filename="POSCAR_ads")

AdsorbateSiteFinder 内部用 Delaunay 三角剖分找对称不等价位点；add_adsorbate 的 height（吸附距离）等参数不同版本略有差异，建议 help(AdsorbateSiteFinder.add_adsorbate) 确认。

7. Materials Project 数据获取

注意：需要 Materials Project API key。 到 materialsproject.org 注册后在 Dashboard 生成。新版 API 与旧 legacy API 是两套体系：

  • 旧：pymatgen.ext.matproj.MPRester + PMG_MAPI_KEY（legacy.materialsproject.org，已弃用）
  • 新（当前推荐）：mp_api.client.MPRester + API key（api.materialsproject.org）

新 API 使用方式：

from mp_api.client import MPRester

API_KEY = "你的_API_KEY"   # 不要硬编码到提交的脚本里，建议放环境变量

with MPRester(api_key=API_KEY) as mpr:
    # 1. 按 material_id 下载结构
    structure = mpr.get_structure_by_material_id("mp-19017")  # 例：LiFePO4 之类

    # 2. 按组分搜索（chemsys 用连字符连接元素）
    docs = mpr.summary.search(
        chemsys="Cu-O",
        fields=["material_id", "formula_pretty",
                "formation_energy_per_atom", "band_gap", "structure"],
    )
    for d in docs[:5]:
        print(d.material_id, d.formula_pretty,
              d.formation_energy_per_atom, d.band_gap)

    # 3. 直接拿相图所需的能量条目
    entries = mpr.get_entries_in_chemsys(["Cu", "O"])
    print("条目数:", len(entries))

常用字段名（mpr.summary.search 的 fields）：material_id、formula_pretty、structure、formation_energy_per_atom、energy_above_hull、band_gap、is_metal、symmetry。

本地可跑替代（无需联网/key）：直接按 §3/§5 手工 Structure.from_spacegroup 构造，或用 §8 自己填能量做相图。日常练习/教学时优先用本地构造，避免网络与 key 干扰。

推荐做法：把 API key 放环境变量，脚本里读取：

import os
from mp_api.client import MPRester

key = os.environ.get("MP_API_KEY")
with MPRester(api_key=key) as mpr:
    ...

8. PhaseDiagram

相图/凸包是判断"某物相是否稳定、形成能多少"的标准工具。核心在 pymatgen.analysis.phase_diagram。

from pymatgen.analysis.phase_diagram import PhaseDiagram, PDPlotter
from pymatgen.entries.computed_entries import ComputedEntry
from pymatgen.core import Composition

# 每个条目 = 某物相 + 其 DFT 总能量（eV）
# 注意：energy 是"整个化学式（原胞）"的总能量，PhaseDiagram 内部会自动换算成每原子能量
entries = [
    ComputedEntry(Composition("Cu"),   -3.7),    # 纯 Cu
    ComputedEntry(Composition("Cu2O"), -11.6),   # Cu2O 原胞（3 原子）
    ComputedEntry(Composition("CuO"),  -7.4),    # CuO 原胞（2 原子）
]

pd = PhaseDiagram(entries)

# 稳定相（落在凸包上的相）
print("=== 稳定相 ===")
for e in pd.stable_entries:
    print(e.composition.reduced_formula, round(e.energy_per_atom, 4))

# 形成能（相对组成元素参考态）
print("=== 形成能 (eV/atom) ===")
for e in entries:
    print(e.composition.reduced_formula, round(pd.get_form_energy_per_atom(e), 4))

# 距离凸包的能差：e_above_hull == 0 表示热力学稳定
print("=== 凸包能差 (eV/atom) ===")
for e in entries:
    decomp, e_above_hull = pd.get_decomp_and_e_above_hull(e)
    print(e.composition.reduced_formula, round(e_above_hull, 4), "->", decomp)

能量凸包稳定性判据：e_above_hull == 0（或 < 某个阈值如 0.05–0.1 eV/atom）视为（亚）稳定；> 0 表示会分解成凸包上相邻相的混合物。get_decomp_and_e_above_hull 返回 (分解反应, 能差)。

绘图（需 matplotlib）：

import matplotlib.pyplot as plt
plotter = PDPlotter(pd)
plotter.get_plot()          # 二元 = 2D 曲线
plt.savefig("pd.png", dpi=150)
# 三元体系会生成 3D 凸包

从 DFT 输出构造条目（比手填能量更真实）：

from pymatgen.entries.computed_entries import ComputedStructureEntry
from pymatgen.io.vasp import Vasprun

entries = []
for name in ["cu", "cu2o", "cuo"]:
    vr = Vasprun(f"{name}/vasprun.xml", parse_potcar_file=False)
    e = ComputedStructureEntry(vr.final_structure, vr.final_energy)
    entries.append(e)

pd = PhaseDiagram(entries)
for e in entries:
    decomp, eh = pd.get_decomp_and_e_above_hull(e)
    print(e.composition.reduced_formula, f"{eh:.3f} eV/atom")

9. 与 VASP 对接

9.1 读 vasprun.xml：能量 / 结构 / DOS / 能带

from pymatgen.io.vasp import Vasprun

vr = Vasprun("vasprun.xml", parse_potcar_file=False)  # 无 POTCAR 时务必关掉，否则会告警/报错

# 能量与结构
E = vr.final_energy                     # 最终总能（eV）
final_struct = vr.final_structure       # 弛豫后结构
print("E =", E, "eV")
print("电子步是否收敛:", vr.converged)
print("离子步是否收敛:", vr.converged_ionic)

# 态密度
dos = vr.complete_dos                   # CompleteDos 对象（含分波态密度）
tdos = vr.tdos                          # 总 DOS
print(tdos.energies.shape, tdos.densities.shape)

# 带隙 / VBM / CBM（单点或非 line-mode 也能给出带隙信息）
gap, cbm, vbm, is_direct = vr.eigenvalue_band_properties
print(f"gap={gap:.3f} eV, VBM={vbm:.3f}, CBM={cbm:.3f}, direct={is_direct}")

能带（需 line-mode KPOINTS 计算）：

bs = vr.get_band_structure(kpoints_filename="KPOINTS", line_mode=True)
print(bs.get_band_gap())       # 带隙
print(bs.bands[Spin.up][0])    # 自旋向上第 0 条能带的能量序列

9.2 读 OUTCAR / OSZICAR

from pymatgen.io.vasp import Outcar, Oszicar

outcar = Outcar("OUTCAR")
print(outcar.final_energy)     # eV
print(outcar.is_stopped)       # 是否正常结束
print(outcar.magnetization)    # 总磁矩（μB）

oszicar = Oszicar("OSZICAR")
print(oszicar.final_energy)
print(oszicar.ionic_steps)     # 每个离子步的能量列表

9.3 生成 VASP 输入

from pymatgen.io.vasp.sets import MPRelaxSet, MPStaticSet

# 用 Materials Project 风格参数生成输入文件（INCAR/POSCAR/POTCAR/KPOINTS）
MPRelaxSet(final_struct).write_input("relax")
MPStaticSet(final_struct).write_input("static")

MPRelaxSet/MPStaticSet 写 POTCAR 需要配置赝势库（见 §2.3）。只做结构后处理则不需要。

9.4 简单 DOS 绘图

import matplotlib.pyplot as plt
from pymatgen.io.vasp import Vasprun

vr = Vasprun("vasprun.xml", parse_potcar_file=False)
tdos = vr.tdos
plt.plot(tdos.energies, tdos.densities[Spin.up], label="up")
plt.plot(tdos.energies, tdos.densities[Spin.down], label="down")
plt.xlabel("E (eV)")
plt.ylabel("DOS")
plt.legend()
plt.show()

Spin 在 pymatgen.electronic_structure.core，常用时 from pymatgen.electronic_structure.core import Spin。

10. 分析工具

10.1 对称性

from pymatgen.symmetry.analyzer import SpacegroupAnalyzer

finder = SpacegroupAnalyzer(s)
print(finder.get_space_group_symbol())      # 空间群符号
print(finder.get_crystal_system())          # 晶系

prim = finder.get_primitive_standard_structure()     # 标准原胞
conv = finder.get_conventional_standard_structure()  # 标准惯用胞

10.2 近邻（配位数、键长）

方法 A：按半径简单判近邻

neighbors = s.get_neighbors(s[0], r=3.0)   # 距 s[0] 3 Å 内的原子
for n in neighbors:
    print(n.species, round(n.nn_distance, 3))   # n.nn_distance 为真实最近邻距离

方法 B：几何算法（更鲁棒，推荐）

from pymatgen.analysis.local_env import CrystalNN

cnn = CrystalNN()
for site in s:
    nn_list = cnn.get_nn_info(s, site.index)
    print(site.specie, "配位数", len(nn_list))
    for nn in nn_list:
        print("   ->", nn["site"].specie, round(nn["length"], 3))   # nn['length'] 为键长

10.3 距离矩阵

dm = s.distance_matrix   # 周期边界下的最短原子间距矩阵 (N,N)
print(dm)

11. 科研实战示例

示例 1：从零构建 Cu(111) 表面并加 CO（CO2 还原催化剂模型）

# build_cu_surface.py
from pymatgen.core import Structure, Lattice
from pymatgen.core.surface import SlabGenerator

# 体相 fcc Cu
cu = Structure.from_spacegroup("Fm-3m", Lattice.cubic(3.615), ["Cu"], [[0, 0, 0]])

# 切 (111) 表面
slabgen = SlabGenerator(cu, (1, 1, 1), min_slab_size=8, min_vacuum_size=15,
                        lll_reduce=True, center_slab=True)
slab = slabgen.get_slab(shift=0)

# 扩成 3x3 表面超胞（催化常用，减小吸附质自相互作用）
slab = slab.make_supercell([[3, 0, 0], [0, 3, 0], [0, 0, 1]])

print("原子数:", slab.num_sites, " 晶格c:", round(slab.lattice.c, 2))
slab.to(filename="POSCAR_slab")
print("已写出 POSCAR_slab")

示例 2：读 VASP 输出算吸附能

吸附能定义（以 CO 吸附为例）：E_ads = E(slab+CO) - E(slab) - E(CO)，负值表示放热、吸附稳定。

# ads_energy.py
from pymatgen.io.vasp import Vasprun

def read_energy(path):
    return Vasprun(path, parse_potcar_file=False).final_energy

E_slab     = read_energy("slab/vasprun.xml")       # 干净表面
E_ads_slab = read_energy("ads_co/vasprun.xml")     # 吸附 CO 后
E_CO       = read_energy("co/vasprun.xml")         # 孤立 CO（放 20 Å 真空盒子算）

E_ads = E_ads_slab - E_slab - E_CO
print(f"CO 吸附能 = {E_ads:.3f} eV")

示例 3：按形成能筛选稳定相（本地 VASP 结果）

# stability.py
from pymatgen.io.vasp import Vasprun
from pymatgen.entries.computed_entries import ComputedStructureEntry
from pymatgen.analysis.phase_diagram import PhaseDiagram

phases = ["cu", "cu2o", "cuo"]   # 各自一个目录，含弛豫后的 vasprun.xml
entries = []
for p in phases:
    vr = Vasprun(f"{p}/vasprun.xml", parse_potcar_file=False)
    entries.append(ComputedStructureEntry(vr.final_structure, vr.final_energy))

pd = PhaseDiagram(entries)
print("稳定相:", [e.composition.reduced_formula for e in pd.stable_entries])
for e in entries:
    decomp, eh = pd.get_decomp_and_e_above_hull(e)
    print(f"{e.composition.reduced_formula:8s} E_hull={eh:6.3f} eV/atom")

示例 4：从 Materials Project 批量下载并切表面（需 API key）

# mp_surface_screen.py  —— 需 Materials Project API key
import os
from mp_api.client import MPRester
from pymatgen.core.surface import SlabGenerator

key = os.environ["MP_API_KEY"]
with MPRester(api_key=key) as mpr:
    # 找 Cu 的稳定氧化物
    docs = mpr.summary.search(
        chemsys="Cu-O",
        energy_above_hull=(0, 0.05),     # 只看稳定/近稳定相
        fields=["material_id", "formula_pretty", "structure"],
    )
    for d in docs[:3]:
        s = d.structure
        # 对每个稳定相切低指数面
        for mi in [(1, 0, 0), (1, 1, 0), (1, 1, 1)]:
            try:
                slab = SlabGenerator(s, mi, min_slab_size=8,
                                     min_vacuum_size=15, center_slab=True).get_slab(shift=0)
                fname = f"{d.material_id}_{''.join(map(str, mi))}.vasp"
                slab.to(filename=fname)
                print("写出", fname)
            except Exception as ex:
                print("跳过", d.material_id, mi, "原因:", ex)

12. 常见坑

  1. Materials Project API key 新旧混淆：pymatgen.ext.matproj.MPRester（legacy）已弃用，务必用 mp_api.client.MPRester（需要 mp-api 包，随新版 pymatgen 安装）。新 key 在 materialsproject.org 的 Dashboard 生成。

  2. Vasprun 报 POTCAR 相关告警/错误：读 vasprun.xml 时若目录无 POTCAR，务必 parse_potcar_file=False。新版该参数默认值有变化，显式写死最稳。

  3. API 版本变动：pymatgen 近年调整过部分接口（如 MPRester 位置、Molecule.substitute/Structure.substitute 签名、Vasprun 参数、SlabGenerator/AdsorbateSiteFinder 细节）。不确定时 help(Class.method) 或查官方文档，勿照搬 2020 年以前的旧教程。

  4. 晶格精度：CIF/POSCAR 互转时小数位数可能导致晶格常数轻微变化。写 CIF 可用 CifWriter(s, significant_figures=8) 提高精度；写 POSCAR 默认足够，但做晶格敏感分析（声子、弹性）要保留充分位数。

  5. Python 3.14 依赖 wheel：部分二进制依赖（spglib/scipy）在新 Python 下可能暂缺预编译包，报编译错误时改 conda-forge 安装或退到 3.12/3.13。

  6. Structure.to 的返回：to(fmt=...) 返回字符串，to(filename=...) 返回 None 并写文件。二者别混用。

  7. slab 真空层：min_vacuum_size 太小会导致相邻周期 slab 相互作用、能量不准；带电/偶极校正体系要更大（≥ 20 Å）并配合 LDIPOL/IDIPOL。

  8. 能量单位：Vasprun.final_energy 是 eV；ComputedEntry 的 energy 是"整个原胞"总能，PhaseDiagram 内部自动按每原子处理，别把 per-atom 值当 total 传。

13. 延伸资源

提示：本文所有代码基于 pymatgen 2025.x 编写，若你本机是 2026.x 或更新版本，个别方法签名可能微调。养成习惯：写脚本前 import pymatgen; print(pymatgen.core.__version__)，对不确定的接口 help() 一眼确认，比任何教程都可靠。