LAMMPS 层间相互作用势 pair_style ilp/graphene/hbn 完全指南：石墨烯/h-BN 异质结 ILP 势的参数、原理与实战配置-开发者社区

科研
科学计算
高性能计算

【免费下载链接】lammps

Public development project of the LAMMPS MD software package

项目地址：https://gitcode.com/gh_mirrors/la/lammps

导读

本文围绕 LAMMPS 中的pair_style ilp/graphene/hbn命令，系统讲解依赖堆叠构型（registry-dependent）的层间相互作用势 ILP（Interlayer Potential）的数学形式、参数文件格式、与hybrid/overlay、coul/shield的联合使用方法，以及通过compute pair分解 vdW 吸引能与排斥能的能量诊断技巧。读完本文，你将能够在 LAMMPS 中完整搭建石墨烯/石墨烯、h-BN/h-BN 及石墨烯/h-BN 异质结的双层乃至多层体系模拟，并理解层法向量计算、taper 函数等底层实现细节。

命令语法与参数

pair_style ilp/graphene/hbn必须在hybrid/overlay中作为子风格使用，语法如下：

pair_style [hybrid/overlay ...] ilp/graphene/hbn cutoff tap_flag

cutoff：全局截断距离（单位与当前units设置一致），示例中为 16.0 Å；
tap_flag：0 或 1，控制是否开启 taper 函数（能量/力的平滑衰减函数）。默认值为1（开启）。

该风格有一个加速变体ilp/graphene/hbn/opt，属于 OPT 包，参数与结果一致（除舍入精度差异外），见文档 accel_styles.rst。加速风格可通过输入脚本中的-suffix opt命令行开关或suffix命令启用。

从源码看，参数解析位于 pair_ilp_graphene_hbn.cpp 的settings()中：参数个数必须为 1 或 2；cut_global直接取自第一个参数；若提供第二个参数则覆盖默认的tap_flag。构造函数 pair_ilp_graphene_hbn.cpp 将tap_flag初始化为 1，并在初始化时检查pair_style必须为^hybrid/overlay，否则报错 "Pair style ilp/graphene/hbn must be used as sub-style with hybrid/overlay"。

势能数学形式与物理原理

核心公式

ILP 势的总能量为：

E = 1/2 * Σ_i Σ_{j≠i} V_ij

其中原子对相互作用：

V_ij = Tap(r_ij) * { exp(-α(r_ij/β - 1)) * [ε + f(ρ_ij) + f(ρ_ji)] - 1/[1 + exp(-d*((r_ij/(s_R·r_eff)) - 1))] * C_6/r_ij^6 }

各项含义：

ρ_ij² = r_ij² - (r_ij · n_i)²、ρ_ji² = r_ij² - (r_ij · n_j)²：横向距离（原子对在各自层平面内的投影距离），正是"依赖堆叠构型"的来源；
f(ρ) = C * exp(-(ρ/δ)²)：层间斥力的横向高斯调制项；
n_i、n_j：原子所在层的法向量，按 Kolmogorov & Crespi 的方法计算（见下文）；
Tap(r_ij) = 20(r_ij/R_cut)⁷ - 70(r_ij/R_cut)⁶ + 84(r_ij/R_cut)⁵ - 35(r_ij/R_cut)⁴ + 1：taper 函数，在截断半径r_c之外提供直至三阶导数连续的平滑截断。

各参数定义见 Leven et al.（J. Chem. Phys. 140, 104106 (2014)）与 Maaravi et al.（J. Phys. Chem. C 121, 22826 (2017)）。

Taper 函数的源码实现

在 interlayer_taper.h 中，calc_Tap与calc_dTap通过霍纳（Horner）形式求值Tap_coeff[8] = {1, 0, 0, 0, -35, 84, -70, 20}对应的七次多项式及其导数；当r_ij/Rcut ≥ 1时 Tap 与 dTap 均为 0。该头文件被 ILP 家族多个势（ilp/graphene/hbn、ilp/tmd、saip/metal、aip/water/2dm）共用，因此 taper 行为在这些层间势中完全一致。

法向量的计算

法向量按 Kolmogorov & Crespi（Phys. Rev. B 71, 235415 (2005)）的方法计算：对每个原子，取其所在层内的近邻原子（同一 molecule id、距离小于参数文件中rcut的原子，见ILP_neigh()中rsq < cutILPsq[...]且atom->molecule[i] == atom->molecule[j]的判定，pair_ilp_graphene_hbn.cpp），构造相邻矢量并求叉积得到法向量：

若近邻数 ≤ 1：法向量直接取(0, 0, 1)；
若近邻数为 2：取两个矢量的叉积；
若近邻数为 3：对三个近邻平面法向量求平均。

同时源码计算了法向量对自身坐标（dnormdri）及各近邻坐标（dnormal）的一阶导数，用于斥力项中涉及法向量的力贡献——这正是该势"多体"特性的体现（manybody_flag = 1，见构造函数）。值得注意的是，若某原子的层内近邻数超过 3，ILP_neigh()会报错 "There are too many neighbors for some atoms, please check your configuration"，因此层内原子数需合理（如六方层内每个原子通常有 3 个近邻）。

关键注意事项

必须包含所有原子对以构建用于计算法向量的邻居列表。由于法向量依赖层内邻居，LAMMPS 会为该风格请求包含 ghost 原子的完整邻居列表（NeighConst::REQ_FULL | NeighConst::REQ_GHOST，见 pair_ilp_graphene_hbn.cpp 的init_style()）。

适用体系与势的"组合使用"要求

该势专用于两层不同的石墨烯、六方氮化硼（h-BN）及其异质结之间的层间相互作用。要获得真实模拟结果，必须与层内（intralayer）势配合使用，例如 AIREBO 或 Tersoff。

为避免影响层内性质，同一层内的原子之间不应有层间相互作用。每个原子需要有"层标识符"：同一层的原子之间用合适的层内势，不同层的原子之间用 ILP。在 LAMMPS 中，molecule id 被选作层标识符，因此使用该势必须使用atom_style full的数据文件（init_style()中会检查atom->molecule_flag否则报错 "requires atom attribute molecule"）。

此外，两个参数文件的差别：

BNCH.ILP（potentials/BNCH.ILP）：Ouyang et al. Nano Lett. 18, 6009-6016 (2018) 发布的新参数，同时适用于短程与长程相互作用，适合高压（层间距小于平衡距离）模拟；
BNCH-old.ILP（potentials/BNCH-old.ILP）：Leven et al.（2016）与 Maaravi et al.（2017）的旧参数，仅适用于长程相互作用区。

参数文件中每一行包含 11 个参数（以 BNCH.ILP 为例）：

# beta(A) alpha delta(A) epsilon(meV) C(meV) d sR reff(A) C6(meV*A^6) S rcut C C 3.205843 7.511126 1.235334 1.528338E-5 37.530428 15.499947 0.7954443 3.681440 25.714535E3 1.0 2.0

其中beta对应公式中的z0，alpha为衰减指数，delta为高斯宽度，epsilon为势阱深度，C为高斯幅度，d控制 vdW 项的开关陡峭度，sR与reff之积构成 vdW 项的特征长度（源码中预计算seff = sR * reff），C6为 vdW 系数，S用于整体缩放能量（源码中meV = 1e-3 * S，S = 43.3634 meV ≈ 1 kcal/mol，读取时对 C、C6、epsilon 乘以此因子，见 pair_ilp_graphene_hbn.cpp），rcut用于构建计算法向量的邻居列表（cutILPsq，见 pair_ilp_graphene_hbn.cpp）。

注意：BNCH.ILP 中 B–N 原子对的参数是按照屏蔽库仑相互作用coul/shield拟合的。因此要正确模拟 h-BN 的性质，该势必须与 pair_style coul/shield 联合使用。

参数文件格式与单位约定

BNCH.ILP 势参数文件按metal 单位拟合，能量单位为meV。文件头包含DATE、UNITS: metal、CONTRIBUTOR与CITATION元信息；每行 11 个数值字段，元素对覆盖 C、B、N、H 之间的全部对称组合（文件同时给出对称部分，源码要求每个元素对必须有且仅有一个条目，重复或缺失都会报错，见 pair_ilp_graphene_hbn.cpp）。

读取逻辑（read_file()）中，只解析pair_coeff中列出的元素组合，其余行跳过。若模拟不使用 metal 单位，需自行准备单位换算后的自定义 BNCH.ILP 文件——参数以 meV 表示，换算时注意能量单位转换（源码通过PotentialFileReader与utils::get_conversion_factor(utils::ENERGY, ...)支持单位自动换算）。

实战配置：石墨烯/h-BN 异质结双层体系

以下配置取自官方示例 in.grhBN（异质结 gr/hBN，碳吸附层与 hBN 膜层）：

# Initialization units metal boundary p p p atom_style full read_data gr_hBN_Cstack_2L_noH.data mass 1 10.8110 # boron mass (g/mole) | membrane mass 2 14.0067 # nitrogen mass (g/mole) | membrane mass 3 12.0107 # carbon mass (g/mole) | adsorbate group hBN molecule 1 group gr molecule 2 pair_style hybrid/overlay rebo tersoff ilp/graphene/hbn 16.0 coul/shield 16.0 pair_coeff * * rebo CH.rebo NULL NULL C # chemical pair_coeff * * tersoff BNC.tersoff B N NULL # chemical pair_coeff * * ilp/graphene/hbn BNCH.ILP B N C # long range pair_coeff 1 1 coul/shield 0.70 pair_coeff 1 2 coul/shield 0.69498201415576216335 pair_coeff 2 2 coul/shield 0.69 neighbor 2.0 bin neigh_modify every 1 neigh_modify delay 0 neigh_modify check yes timestep 0.001 velocity all create 300.0 12345 dist gaussian mom yes rot yes fix thermostat all nve compute 0 all pair rebo compute 1 all pair tersoff compute 2 all pair ilp/graphene/hbn compute 3 all pair coul/shield variable REBO equal c_0 variable Tersoff equal c_1 variable EILP equal c_2 # total interlayer energy variable Evdw equal c_2[1] # attractive energy variable Erep equal c_2[2] # repulsive energy variable Ecoul equal c_3 thermo 100 thermo_style custom step etotal pe ke v_Evdw v_Tersoff v_REBO v_EILP v_Erep v_Ecoul temp run 1000

要点解读：

hybrid/overlay 同时叠加三个子势：rebo（石墨烯层内 C–C 化学键）、tersoff（hBN 层内 B–N 化学键）、ilp/graphene/hbn（跨层长程相互作用）、coul/shield（B–N 屏蔽库仑）。pair_coeff中用NULL关闭不需要的元素组合，使各子势"各管一摊"；
不同层用不同 molecule id（group hBN molecule 1、group gr molecule 2），ILP 只作用于molecule[i] != molecule[j]的原子对（pair_ilp_graphene_hbn.cpp 与 L544 的判断）；
双层石墨烯示例 in.bilayer-graphene 使用hybrid/overlay rebo ilp/graphene/hbn 16.0且pair_coeff * * ilp/graphene/hbn BNCH.ILP C C；双层 hBN 示例 in.bilayer-hBN 使用tersoff + ilp/graphene/hbn + coul/shield组合，其中 coul/shield 的屏蔽长度取 0.70（B–B）、0.69498（B–N）、0.69（N–N）；
完整示例目录见 examples/PACKAGES/interlayer/ilp_graphene_hbn，内含 4 个输入脚本、4 个 data 文件、势文件及运行日志，可直接参考或对照验证。

能量分解与热力学输出

该 pair style 会把总层间势能分解为两个子类别，通过 compute pair 以长度为 2 的向量形式访问：

E_vdW：vdW（吸引）能量；
E_Rep：排斥能量。

在输入脚本中加入如下命令即可将两项写入 log 文件（带列标题）：

compute 0 all pair ilp/graphene/hbn variable Evdw equal c_0[1] variable Erep equal c_0[2] thermo_style custom step temp epair v_Erep v_Evdw

源码中pvector长度为 2（构造函数nextra = 2），calc_FvdW()累加pvector[0]（吸引项Vilp = -C6*r6inv/TSvdw，pair_ilp_graphene_hbn.cpp），calc_FRep()累加pvector[1]（排斥项Vilp = exp(-λ(r-z0)) * (0.5ε + f(ρ))，pair_ilp_graphene_hbn.cpp），两者均乘上 taper 函数后分别统计。

兼容性、限制与注意事项

支持的选项

该风格不支持pair_modify的 mix、shift、table、tail 选项；
该风格不将参数信息写入二进制 restart 文件（restartinfo = 0），因为参数存于势文件中；读取 restart 文件后必须重新指定pair_style与pair_coeff。

构建要求

属于INTERLAYER 包，仅当 LAMMPS 以该包构建时启用。构建方式参见 Build package 文档，例如-D PKG_INTERLAYER=yes；OPT 加速变体需要同时启用 OPT 包；
该风格要求newton对 pair 相互作用为 on（默认即为 on），否则init_style()报错 "requires newton pair on"；
必须使用atom_style full（需要 molecule 属性作为层标识）；
必须在hybrid/overlay中作为子风格使用；其他不需要的相互作用可用pair_style none置零。

测试验证

单元测试 manybody-pair-ilp-graphene-hbn.yaml（tap_flag=1，48 原子双层体系，初始 vdW 能量 -1.3859）与 manybody-pair-ilp-graphene-hbn_notaper.yaml（tap_flag=0，需comm_modify cutoff 16.0，初始 vdW 能量 -1.7414）分别验证了 taper 开启/关闭两种模式下的能量、应力与力，可作为正确性参照。注意notaper变体要求通信截断不小于全局截断，否则远距离原子对无法正确参与相互作用。