LAMMPS计算气体分子在固体表面吸附能的脚本编写

LAMMPS吸附能计算的基本思路

与VASP通过量子力学波函数计算能量不同，LAMMPS用经典力场计算原子间相互作用能。吸附能的公式相同：

E_ads = E(total) – E(surface) – E(adsorbate)

但实现方式完全不同——LAMMPS中需要编写in脚本，通过group定义和compute命令提取各部分能量。关键挑战在于力场中很多能量项（如pair_energy）是全局计算的，无法直接分配到特定原子组上。本项目在实践中总结了一套可靠的能量分离方法。更多[LAMMPS相关计算](https://www.keyanxueshu.com/category/md/)的案例可参考站内文章。

能量分离的三种方法

方法1：三次独立计算

最直观的方法：建立三个独立的data文件，分别跑三个计算。

优点：简单可靠，不存在能量分配问题。

缺点：需要准备三个data文件，计算B和C的构型必须与A中完全一致（不额外优化）。

本项目的标准做法：计算A完成后，从最终构型中”删除”吸附质原子得到计算B的初始构型，”删除”表面原子得到计算C的初始构型。不重新优化，直接做单点能量计算。

方法2：compute group/group

LAMMPS内置的compute group/group命令可以直接计算两个group之间的相互作用能：

compute 1 surface/group adsorbate group/group

这给出表面和吸附质之间的相互作用能（不含各自内部能量），就是吸附能的负值（符号相反）。

关键限制：compute group/group只计算pair_style中的短程相互作用。如果力场包含长程项（如KSPACE Ewald），需要加kspace yes关键词。但即使如此，KSPACE中的自相互作用校正可能不完整。

方法3：能量分解（pair epair）

通过将不同能量项分配到不同group：

compute 1 adsorbate pe/atom pair
compute 2 surface pe/atom pair

理论上，总pair energy = surface内部 + adsorbate内部 + surface-adsorbate之间。前两项可以通过compute获取，后者用差值。但这个方法的精度取决于力场实现是否正确支持group级别的能量分解。

本项目推荐方法1（三次独立计算）作为主方法，方法2（compute group/group）作为交叉验证。

完整in脚本实例

以N₂分子在石墨烯表面的物理吸附为例：

计算A：表面+吸附质

# ———- 基本设置 ———-
units           metal
atom_style      full
boundary        p p p

# ———- 读取结构 ———-
read_data       graphene_n2.data

# ———- 力场设置 ———-
pair_style      lj/cut 10.0
pair_coeff      1 1 0.0029 3.28   # C-C (石墨烯)
pair_coeff      2 2 0.0066 3.31   # N-N (N₂)
pair_coeff      1 2 0.0044 3.30   # C-N (混合)

# ———- 定义组 ———-
group           carbon type 1
group           nitrogen type 2
group           n2_molecule id 1201 1202   # N₂的两个原子
group           surface id 1:1200           # 石墨烯碳原子
group           all_atoms type 1 2

# ———- 最小化 ———-
neighbor        2.0 bin
neigh_modify    every 1 delay 0 check yes

# 固定底层石墨烯
fix             1 surface setforce 0.0 0.0 0.0

# N₂分子键约束
fix             2 n2_molecule rigid/npt small molecule

# 最小化
min_style       cg
minimize        1.0e-8 1.0e-8 1000 10000

# ———- 输出能量 ———-
variable        E_total equal etotal
variable        E_pair equal epair
variable        temp equal temp

thermo          100
thermo_style    custom step v_E_total v_E_pair v_temp pe ke etotal press

# 交互作用能（交叉验证）
compute         int_energy surface group/group n2_molecule
variable        E_interaction equal c_int_energy

# ———- 运行 ———-
dump            1 all custom 1000 dumpA.lammpstrj id type x y z
run             0    # 单点能量计算

print           “E_total = ${E_total}”
print           “E_interaction = ${E_interaction}”

计算B：纯表面

从计算A的最终构型中删除N₂原子，生成graphene_only.data，运行：

# 同样的力场设置和最小化
# 只输出总能量
run 0
print           “E_surface = ${E_total}”

计算C：孤立N₂

构建一个含2个N原子的盒子（≥15Å边长）：

units           metal
atom_style      atomic
boundary        p p p
read_data       n2_isolated.data

pair_style      lj/cut 10.0
pair_coeff      1 1 0.0066 3.31

neighbor        2.0 bin
min_style       cg
minimize        1.0e-10 1.0e-10 1000 10000

run 0
print           “E_adsorbate = ${E_total}”

吸附能计算

E_ads = E_A – E_B – E_C

本项目在N₂-graphene体系中的计算结果：

与DFT-D3计算值-0.135 eV的偏差为9%，对于经典力场来说是可以接受的精度。

力场选择指南

LAMMPS吸附能计算的精度高度依赖力场质量：

混合力场处理：当表面和吸附质属于不同力场体系时（如金属表面+有机分子），使用hybrid或hybrid/overlay：

pair_style      hybrid eam/alloy lj/cut 10.0
pair_coeff      1 1 eam/alloy Cu_u3.eam.alloy   # Cu-Cu (EAM)
pair_coeff      2 2 lj/cut 0.0066 3.31           # N-N (LJ)
pair_coeff      1 2 lj/cut 0.0044 3.30           # Cu-N (LJ混合规则)

混合规则（Lorentz-Berthelot）：σ_12 = (σ_1 + σ_2)/2, ε_12 = √(ε_1 × ε_2)

但混合规则不总是可靠，特别是金属-有机界面。本项目建议从DFT计算中拟合交叉参数，而非直接用混合规则估算。关于[科研计算服务](https://www.keyanxueshu.com/)中力场拟合的方法可参考站内文章。

周期性边界与截断半径

截断半径的陷阱：LJ势在截断处不连续，会导致能量跳跃。解决方案：

1. 使用pair_modify shift yes——将势能曲线平移到截断处为零（能量连续但力不连续）
2. 使用pair_style lj/smooth——在截断附近做平滑过渡（力和能量都连续）
3. 使用pair_style lj/cut/coul/long配合kspace——适合含电荷体系

温度效应：从0K到有限温度

上述脚本计算的是0K下的静态吸附能。如果要考虑温度效应：

# NVT系综下运行MD
fix             3 all_atoms nvt temp 300.0 300.0 0.1

# 每N步计算一次交互作用能
compute         int_energy surface group/group n2_molecule
fix             4 all_atoms ave/time 100 10 1000 v_E_interaction file interaction.dat

run             100000    # 100 ps

然后对interaction.dat中的数据取平均和方差，得到有限温度下的平均吸附能及其热涨落。本项目发现，300K下物理吸附能的热涨落可达0.05-0.10 eV，约占吸附能本身的30-50%。

常见错误

LAMMPS计算吸附能的优势在于速度快、体系规模大（可处理数万原子），适合做参数扫描和温度效应分析。但精度受限于力场质量——对于化学吸附或涉及电荷转移的体系，仍然需要DFT或AIMD来获取可靠结果。本项目推荐的方法是多级策略：LAMMPS粗筛→DFT精算，兼顾效率和精度。