MS表面能计算：从建模到结果的实战经验

表面能是材料表面科学中最基础也最重要的物理量之一。无论是催化、薄膜生长还是润湿性研究，表面能计算都是第一步。用Materials Studio的CASTEP模块算表面能看似简单——建个slab模型跑个优化就行了——但实际操作中有很多细节决定结果是否可靠。本文分享我在大量表面能计算项目中积累的经验。

一、表面能的物理意义与公式

1.1 基本公式

表面能γ的定义：

γ = (E_slab – N × E_bulk) / (2 × A)

其中：

• E_slab：slab模型的总能量
• E_bulk：体相单原子能量
• N：slab中的原子数
• A：一个表面的面积
• 因子2：slab有上下两个表面

1.2 不同晶面的表面能

以面心立方(FCC)金属为例，常见低指数表面的表面能排序：

γ(111) < γ(100) < γ(110)

这与密排程度一致——密排面断键最少，表面能最低。

经验参考值（PBE计算，单位mJ/m²）：

注意：DFT计算值通常比实验值低20-30%，这是因为实验中表面有重构、污染等因素。

二、MS建模实操经验

2.1 slab厚度选择

这是表面能计算中最关键的一个参数。

经验法则：

• 金属：至少4层，底层2层固定
• 氧化物：至少6层，底层3层固定
• 二维材料：单层即可，无需slab

收敛测试方法：

1. 从3层开始，逐步增加到4、5、6、7层
2. 计算每个厚度下的表面能
3. 当相邻厚度间表面能变化<10 mJ/m²时认为收敛

踩坑经验：我早期算RuO2(110)表面能时用了4层slab，结果比文献值高了15%。后来加到6层才收敛。氧化物的层间相互作用比金属强，需要更厚的slab。

2.2 真空层设置

• 最小值：15 Å（有偶极校正）
• 推荐值：20 Å（无偶极校正时更安全）
• 偶极校正：对极性表面（如氧化物(0001)面）必须开启

在MS中设置：CASTEP → Calculation → Electronic → 勾选”Dipole correction”

2.3 超胞大小

表面能计算不需要大超胞，(1×1)即可：

• 超胞大小只影响k点选择，不直接影响表面能
• 但如果要算吸附，需要(2×2)或(3×3)
• 表面能计算用(1×1)+高密k点，效率最高

2.4 底层原子固定

在MS中：

1. 选中底层2层原子
2. Properties → Constraints → 勾选”Fix Cartesian position”
3. 只优化顶层原子

经验：固定底层模拟体相环境，这是slab模型的基本假设。如果底层也放松，表面能会偏低。

三、计算参数设置经验

3.1 CASTEP参数推荐

# 关键参数设置
Task: Geometry Optimization
Quality: Fine
Energy cutoff: 500 eV（金属400够，氧化物500+）
k-point spacing: 0.03 1/Å
SCF tolerance: 1.0E-6 eV/atom
Max force: 0.03 eV/Å
Max displacement: 0.001 Å
Max iterations: 200

# 电子展宽
Smearing: Gaussian
Smearing width: 0.1 eV（金属）

3.2 k点选择

slab模型沿真空层方向用1个k点：

验证：k点加倍后表面能变化<5 mJ/m²即收敛。

四、完整计算流程

Step 1: 体相优化

先优化体相结构，获得E_bulk：

1. 导入晶体结构（从MS结构库或CIF文件）
2. CASTEP几何优化（ISIF=3，即同时优化晶胞和原子）
3. 记录优化后的总能量E_bulk_total和原子数N_bulk
4. 计算E_bulk = E_bulk_total / N_bulk

关键：体相优化必须用与slab计算完全一致的参数（截断能、泛函、k点密度），否则引入系统误差。

Step 2: 构建slab模型

1. Build → Surfaces → Cleave Surface
2. 选择晶面（如(111)）
3. 设置slab厚度（Fractional depth）
4. 添加真空层（Vacuum = 20 Å）
5. Build → Symmetry → Supercell（如需超胞）

Step 3: slab优化

1. 固定底层原子
2. CASTEP几何优化（只优化原子位置，不优化晶胞）
3. 记录优化后的总能量E_slab

Step 4: 计算表面能

γ = (E_slab – N_slab × E_bulk) / (2 × A)

其中A是slab面内的表面积，可以从MS的Lattice参数算出。

五、特殊情况处理

5.1 非对称slab

有些表面必须用非对称slab（如极性氧化物表面），此时上下表面不同，公式变为：

γ = (E_slab – N × E_bulk – γ_bottom × A) / A

需要先用对称slab计算底面表面能，再代入算顶面。

5.2 表面重构

某些表面在优化过程中会发生重构（如Si(111)-(7×7)），这时：

• 需要更大的超胞
• 优化步数要增加（NSW=500+）
• POTIM要减小（0.2-0.3）
• 结果可能与(1×1)理想表面差别很大

5.3 自旋极化

对磁性材料（Fe、Co、Ni、Cr等）：

• 必须开启自旋极化
• 初始磁矩设置：体相值×1.2（表面增强效应）
• 表面层的磁矩通常比体相大20-50%

六、常见问题与排查

Q1：表面能为负值？

原因：slab不够厚或体相能量不准 解决：增加slab层数，确认体相优化参数一致

Q2：表面能比文献值高很多？

可能原因：

1. slab太薄（<4层）
2. 底层没有固定
3. k点不够密
4. 截断能太低
5. 没有色散校正（DFT-D3）

Q3：优化不收敛？

尝试：

1. 减小POTIM到0.2
2. 增加NSW到500
3. 检查初始结构是否合理
4. 尝试不同的优化算法（CG→BFGS）

Q4：不同晶面表面能差异不显著？

正常情况：对于某些立方晶体，(111)和(100)差异本就不大（<5%）。如果差异<1%，检查是否结构对称性导致。

七、项目报价参考

结语

表面能计算看似简单，但要做好需要严格收敛测试和细致的参数控制。建议新手先用Pt(111)练手，与文献值对比验证后再做自己的体系。如有计算需求或技术问题，欢迎联系我们的计算团队。