pymatgen Wulff构型实操

下载主要模块实操包，解压后在包的根目录执行：

bash run.sh
ls analysis
ls figures

脚本使用集群已有的pymatgen 2026.9.24环境，包含绘图依赖，单次使用1个CPU核心；不启动新的VASP计算。也可以分步执行：

module load pymatgen/2026.9.24
python modules.py
python plot.py

1. 先有表面能，再做形状构造#

WulffShape接收晶格、Miller指数和对应表面能。本例从已有Cu体相与对称Cu(111)薄层的原始输出重新提取数据：

两面等价的对称薄层使用：

gamma = (E_slab - N_slab*E_bulk/N_bulk) / (2*A)
1 eV/Ų = 16.021766208 J/m²

面积从原始CONTCAR晶格矢量叉乘得到，能量使用同轮XML的E0。体相与薄层原始INCAR、POSCAR、KPOINTS、CONTCAR和输出都在附件中；这一步不重新提交任务。该薄层只是教学示例，正式表面能仍需测试厚度、真空和k点等条件。

cat inputs/cu-surface-symmetric/INCAR
cat inputs/cu-surface-symmetric/KPOINTS
cat analysis/RESULTS.json

2. 本例只提供{111}面族#

shape = WulffShape(cu.lattice, [(1, 1, 1)], [gamma])

立方对称性生成等价{111}晶面，得到如下受限构型：

这里没有输入{100}、{110}和其他可能暴露面的能量，因此不能把八面体图称为已经预测了完整的Cu平衡形貌。图中坐标经过归一化，也不能直接读作纳米颗粒尺寸。

3. 用于自己的材料#

先确定相与晶格，再对准备比较的晶面分别完成表面模型及一致参照的能量计算。把各自Miller指数和表面能成对传给WulffShape，注意J/m²与eV/Ų不要混用。

本例只有一个对称面族，面积占比为1；这来自输入集合的限制，不是经过全部候选面的比较后证明其他面不出现。实际吸附环境、化学势或不同终止面也可能改变表面自由能，不能由这里的静态清洁表面能直接推断。pymatgen Wulff接口。