LAMMPS计算吸附能：我的完整实操流程与关键参数

吸附能是衡量吸附剂与吸附质相互作用强度的核心指标，结合自己做过的多个表面吸附体系，聊聊LAMMPS计算吸附能的完整流程与常见问题。

【1. 吸附能的定义：搞清楚在算什么】

吸附能(E_ads)的定义是：E_ads = E_total – (E_surface + E_adsorbate)，即吸附体系的总能量减去孤立基底表面能量和孤立吸附质能量。负值代表放热吸附，绝对值越大，吸附越强。这个公式看起来简单，但实际操作中，每一项的计算条件必须严格一致——基底晶格常数、吸附质构型、边界条件，任何不一致都会引入误差。

我在做第一个气体吸附项目时，基底和吸附质的能量是用不同构型分别算的，后来发现吸附质的初始构型差了几度，引入了几十meV的系统误差。这让我后来养成习惯：所有能量计算用同一套赝势和参数，基底构型从同一个优化好的结构里切。

【2. 能量计算方法】

LAMMPS计算吸附能主流有三种方法：

**直接法(Energy minimization)**：最简单直接，用fix minimize命令分别计算三种构型的能量。做静态计算，不涉及温度和动力学，适合研究平衡吸附构型。我大部分时候用这个方法。

**分子动力学+热力学积分**：考虑温度和振动熵的贡献，适合研究高温吸附行为。计算量大，但能给出更接近实验条件的结果。我在金属表面氢脆研究里用过一次，高温下吸附能比0K直接法低了约15%。

**Metadynamics增强采样**：处理强吸附或弱吸附的极端情况，标准MD可能无法充分采样的体系。这个方法我没在吸附能计算里用过，但在研究表面扩散的时候用过。

【3. 我的标准流程：表面建模 → 能量计算】

第一步是切表面：用LAMMPS的delete_atoms命令或者VESTA/MS从块体结构切出指定晶面(我一般用最密排面，比如FCC的(111)面或者Cu的(100)面)。表面需要足够大(至少5×5×n个原子)，避免吸附质与镜像相互作用。

第二步是固定基底：只优化表面附近的原子，底层原子固定，防止表面弛豫把真实的几何结构搞乱。

“`

# LAMMPS 能量最小化示例

units metal

boundary p p f

atom_style atomic

# 读取结构

read_data surface.data

# 势函数(EAM for metals)

pair_style eam/alloy

pair_coeff * * Cu_zhou.eam.alloy Cu

# 固定底层原子(ID 1-500)

group fixed ID <= 500

fix 1 fixed setforce 0 0 0

# 能量最小化

minimize 1.0e-6 1.0e-8 10000 100000

“`

【4. 几种常见体系的经验】

**金属表面气体吸附**(H₂、O₂、CO)：我一般用EAM势描述金属基底，吸附质用LJ势或MEAM势描述。这类体系的吸附能通常在-0.5 ~ -2.0 eV范围，和DFT结果对比时需要做basis set superposition error(BSSE)校正。

**氧化物表面吸附**：这类体系最难，氧化物表面有极性终止和化学终止的差异，表面模型本身就有多种可能。我在处理氧化物的时候，通常参考文献里已有的表面模型，而不是自己从零构建。

**二维材料(Graphene、MoS₂)吸附**：我用过reaxFF和LJ势两种方案。reaxFF精度高但参数难找，LJ势简单但精度有限。Graphene上的气体吸附(CO₂、CH₄等)用LJ势就够了，吸附能很弱(几十meV)，体系本身是范德华相互作用主导。

【5. 常见问题与误差来源】

表面盒子的尺寸不够是新手常犯的错误：吸附质和它的周期性镜像相互作用，导致结果严重偏离。解决办法是确保盒子在z方向至少有30 Å的真空层。另外，势函数的选取对吸附能的影响很大——同一基底用不同势函数，吸附能可能差出30%以上。重要结论一定要用两种以上的势函数做交叉验证。