☰
LAMMPS 角度势函数 angle_style class2 系列:COMPASS/ClassII 力场交叉项与键解离建模全指南
2026/10/15 15:39:11 网站建设 项目流程
  • 科研
  • 科学计算
  • 高性能计算

【免费下载链接】lammps

Public development project of the LAMMPS MD software package

项目地址:https://gitcode.com/gh_mirrors/la/lammps
点击查看免费下载

导读

angle_style class2是 LAMMPS 中实现 COMPASS(Class II)力场角度项的官方风格,它除了标准的四次多项式角度势外,还引入了 bond-bond 与 bond-angle 两类交叉项,用于准确描述分子内振动的非谐性与键长-键角耦合。本文基于官方文档 doc/src/angle_class2.rst,系统讲解class2、class2xe(支持键解离的指数型扩展)、class2/p6(六阶展开)三种风格的能量公式、系数含义、输入脚本与 data 文件配置方法,并结合 CLASS2/MOFFF 包源码与单元测试给出源码级佐证。读完本文,你将能独立完成 COMPASS/PCFF 类力场角度项的建模,并理解如何用class2xe模拟材料断裂过程。


一、三个角风格速览与适用场景

LAMMPS 在同一系列下提供了三个递进的角风格:

风格全称含义所属包用途
class2COMPASS Class II 标准形式CLASS2固定键长体系的标准 Class II 角度项
class2xeClass II + cross-term + exponentialCLASS2允许键拉伸解离的指数型交叉项,配合 Morse 键势
class2/p6class2 角度项扩展到六阶MOFFF需要更高阶非谐角度展开的体系

三个风格共享同一组 bond-bond 与 bond-angle 交叉项骨架,区别在于角度主项(E_a)的次数以及交叉项是二次型还是 Morse 指数型。其中class2与class2xe需要 LAMMPS 以 CLASS2 包构建,class2/p6需要 MOFFF 包(见 Build_package)。

加速变体:class2/kk、class2xe/kk、class2/p6/kk(KOKKOS 包)以及class2/omp、class2xe/omp、class2/p6/omp(OPENMP 包)均提供与串行版本一致的系数约定,可直接用于多核/GPU 加速运行。


二、class2:COMPASS 角度势的标准形式

2.1 能量公式

class2角度风格采用三项求和:

$$E = E_a + E_{bb} + E_{ba}$$

其中角度主项为四次多项式:

$$E_a = K_2(\theta-\theta_0)^2 + K_3(\theta-\theta_0)^3 + K_4(\theta-\theta_0)^4$$

交叉项分别为:

$$E_{bb} = M (r_{ij}-r_1)(r_{jk}-r_2)$$

$$E_{ba} = N_1(r_{ij}-r_1)(\theta-\theta_0) + N_2(r_{jk}-r_2)(\theta-\theta_0)$$

这里θ₀是平衡角度,r₁、r₂分别是 i–j 与 j–k 两条键的平衡键长。E_bb描述两条键长的耦合,E_ba描述键长与键角的耦合,二者共同补偿单纯角度势无法表达的分子内振动耦合效应。该形式源自 COMPASS(Class II)力场,见 Sun, J Phys Chem B 102, 7338-7364 (1998)(文档引用见 doc/src/angle_class2.rst 末尾(Sun)条目)。

2.2 四类系数的完整定义

class2共需要三类系数设置,分别通过angle_coeff命令的不同分支给出。

①E_a公式的 4 个系数(默认分支,无前缀关键字):

系数单位说明
θ₀度(degrees)平衡角度,LAMMPS 内部自动转换为弧度
K₂能量二次项系数,等效为 energy/radian²
K₃能量三次项系数,等效为 energy/radian³
K₄能量四次项系数,等效为 energy/radian⁴

②E_bb公式的 4 个系数(bb前缀):

系数单位说明
bb—标记为 BondBond 系数的关键字
M能量/距离²键-键耦合强度
r₁距离i–j 键平衡键长
r₂距离j–k 键平衡键长

③E_ba公式的 5 个系数(ba前缀):

系数单位说明
ba—标记为 BondAngle 系数的关键字
N₁能量/距离i–j 键与角度耦合强度
N₂能量/距离j–k 键与角度耦合强度
r₁距离i–j 键平衡键长
r₂距离j–k 键平衡键长

注意:E_ba公式中不单独指定θ₀,它复用E_a公式中的同一平衡角度。

2.3 输入脚本写法

angle_style class2 angle_coeff * 75.0 25.0 0.3 0.002 # 所有角度类型的 E_a 系数 angle_coeff 1 bb 10.5872 1.0119 1.5228 # 类型 1 的 BondBond 系数 angle_coeff * ba 3.6551 24.895 1.0119 1.5228 # 所有类型的 BondAngle 系数

关键细节:bb分支在输入脚本中列出 4 个系数(bb+ M + r₁ + r₂),ba分支列出 5 个系数(ba+ N₁ + N₂ + r₁ + r₂)。但在 data 文件中必须去掉bb/ba关键字:

  • BondBond 系数放在BondBond Coeffs标题下,每行仅列角度类型 + 3 个系数;
  • BondAngle 系数放在BondAngle Coeffs标题下,每行仅列角度类型 + 4 个系数。

这与 angle_coeff 文档中特别注明的“class2 是例外,输入脚本需额外参数指定交叉项系数”的约定一致。

2.4 源码级验证

从 src/CLASS2/angle_class2.cpp 可以看到系数解析逻辑与公式完全对应:

  • coeff()中根据第二个参数是否为"bb"或"ba"分派三类系数存储,bb分支校验narg != 5、ba分支校验narg != 6(即含类型号),其余默认分支校验narg != 5(含类型号 + 4 个系数),与文档表格逐位对应;
  • 角度主项能量按k2*dtheta2 + k3*dtheta3 + k4*dtheta4累加,交叉项按bb_k*dr1*dr2与ba_k1*dr1*dtheta + ba_k2*dr2*dtheta累加(angle_class2.cpp 中compute()与single());
  • write_data()将θ₀以theta0/MY_PI*180.0(弧度转度)写出,并分别输出BondBond Coeffs(3 个系数)与BondAngle Coeffs(4 个系数)节,印证了 data 文件格式;
  • extract()支持按k2/k3/k4/theta0关键字提取内部数组,供compute等模块复用。

单元测试 unittest/force-styles/tests/angle-class2.yaml 展示了完整的多类型系数设置示例:

angle_coeff: | 1 110.1 75.0 5.2 1.0 2 111.0 45.0 2.0 0.5 * bb 0.0 1.0 1.0 2 bb 10.0 1.5 1.5 * ba 0.0 0.0 1.5 1.5 1 ba 20.0 0.0 1.5 1.5

可以看到*通配符可为所有角度类型设置同一组系数,再对个别类型单独覆盖——这正是 LAMMPS 输入脚本中通配符 + 覆盖的推荐用法。

2.5 I/J/K 原子顺序的重要性

一个容易被忽视但至关重要的约束:data 文件Angles节中每个角度列表的 I、J、K 原子顺序,必须与r₁、r₂系数的定义顺序一致。因为E_bb与E_ba公式中的r_ij由 I、J 原子计算,r_jk由 J、K 原子计算。若顺序颠倒,交叉项会作用于错误的键,导致势能面错误。源码中compute()正是按delx1 = x[i1]-x[i2]计算第一条键、delx2 = x[i3]-x[i2]计算第二条键,顺序一旦错位交叉项便张冠李戴。


三、class2xe:引入键解离的指数型交叉项

3.1 为什么需要 class2xe

COMPASS/PCFF 等固定键力场的四次键势与谐波交叉项在键长大幅拉伸时物理性不足,无法描述键解离。ClassII-xe 形式(Kemppainen, npj Computational Materials 11, 341 (2025),见 Howto_bioFF)借鉴 Morse 势把“谐波键 → 可解离键”的改造思路,同步应用于交叉项:x 指 cross-term(交叉项),e 指 exponential(指数函数)。这样键拉伸、键角、二面角、异常角的耦合势都允许键解离,从而支撑材料拉伸、屈服直至断裂的模拟。

3.2 能量公式

E_a = K_2(\theta-\theta_0)^2 + K_3(\theta-\theta_0)^3 + K_4(\theta-\theta_0)^4

角度主项不变,交叉项换成 Morse 型:

$$E_{bb} = D \left[ 1 - e^{-\alpha (r_{ij}-r_1)} \right] \left[ 1 - e^{-\alpha (r_{jk}-r_2)} \right]$$

$$E_{ba} = D_1 \left[ 1 - e^{-\alpha_1 (r_{ij}-r_1)} \right] (\theta-\theta_0) + D_2 \left[ 1 - e^{-\alpha_2 (r_{jk}-r_2)} \right] (\theta-\theta_0)$$

D、D₁、D₂为解离能(dissociation energy),α、α₁、α₂为控制指数衰减速率的常数(量纲为 1/距离)。当键长远大于平衡键长时,(1-e^{-αΔr})趋于 1,交叉项强度趋近常数,实现键解离后的自然衰减。

3.3 系数定义

E_a的 4 个系数与 class2 完全相同(θ₀、K₂、K₃、K₄),θ₀同样以度输入、内部转弧度。

E_bb的 5 个系数(bb前缀,输入脚本每行 5 个):

系数单位说明
bb—关键字
D能量键-键交叉项解离能
α1/距离指数衰减常数
r₁距离i–j 键平衡键长
r₂距离j–k 键平衡键长

data 文件中BondBond Coeffs节每行 = 角度类型 + 4 个系数(去掉bb)。

E_ba的 7 个系数(ba前缀,输入脚本每行 7 个):

系数单位说明
ba—关键字
D₁能量i–j 键-角交叉项解离能
D₂能量j–k 键-角交叉项解离能
α₁1/距离i–j 侧指数衰减常数
α₂1/距离j–k 侧指数衰减常数
r₁距离i–j 键平衡键长
r₂距离j–k 键平衡键长

data 文件中BondAngle Coeffs节每行 = 角度类型 + 6 个系数(去掉ba)。

源码 src/CLASS2/angle_class2xe.cpp 中coeff()的分支校验(bb分支narg != 6、ba分支narg != 8,即含类型号)、compute()中的exp(-bb_alpha*dr)能量累加、write_data()输出的 4 个/6 个系数的 data 节,均与上述约定一一吻合。单元测试 unittest/force-styles/tests/angle-class2xe.yaml 给出了真实感的多类型示例:

angle_coeff: | 1 110.1 41.453 -10.604 5.129 * bb 90.4638 0.1368258 1.53 1.101 2 bb 116.34805 0.22772772 1.501 1.417 * ba 82.6961 98.2315 0.35423657 0.24110803 1.53 1.101 1 ba 82.6961 150.0 0.53340706 0.32185452 1.501 1.417

3.4 断裂后的并行通信限制与补救

重要限制:class2xe只描述键拉伸的解离,一旦键解离并被拉伸到超过并行处理器通信截断距离(processor communication cutoff),并行模拟会因 “atoms missing” 错误崩溃。这通常发生在材料断裂之后,因此要模拟断裂后现象,必须先移除已解离键及与之关联的高阶相互作用(角度、二面角、异常角)。可用以下命令配合:

  • fix bond/react:可在键断裂时重组/移除相互作用,支持通过stabilization关键字对断裂点局部原子使用 fix nve/limit 缓和能量不连续;
  • fix bond/break:直接断开键并删除相关高阶项,但注意其引入的能量景观不连续可能造成不稳定。

按 Howto_bioFF 的说明,PCFF-xe 在不使用上述 fix 时可模拟材料直到断裂,用于计算拉伸/剪切模量、屈服强度、泊松比等力学性质;断裂后现象则必须依赖这两个 fix。

3.5 完整的 class2xe 运行示例

仓库中提供了可直接运行的真实示例 examples/PACKAGES/class2xe/in.NVT_class2xe.script,其完整配置展示了 class2xe 与整个 ClassII-xe 力场的配套使用:

units real atom_style full bond_style morse angle_style class2xe dihedral_style class2xe improper_style class2 special_bonds lj/coul 0 0 1 kspace_style pppm 1.0e-4 pair_style lj/class2/coul/long 12.0 pair_modify mix sixthpower neighbor 2.0 bin neigh_modify delay 0 every 1 check yes one 5000 page 100000 read_data detda_PCFF-2xe.data

该脚本配套的 data 文件 examples/PACKAGES/class2xe/detda_PCFF-2xe.data 由 LUNAR 工具构建(README 见 examples/PACKAGES/class2xe/README),演示了从 PCFF 母力场转换到 PCFF-xe 的完整链路:LUNAR 的atom_typing/all2lmp模块构建 COMPASS 或 PCFF 模型,auto_morse_bond_update模块将其重参数化为 COMPASS-xe/PCFF-xe。注意special_bonds lj/coul 0 0 1是该功能形式的标准配置。


四、class2/p6:角度主项六阶展开

4.1 能量公式

class2/p6将 class2 的角度主项从四次多项式扩展到六次:

$$E_a = K_2(\theta-\theta_0)^2 + K_3(\theta-\theta_0)^3 + K_4(\theta-\theta_0)^4 + K_5(\theta-\theta_0)^5 + K_6(\theta-\theta_0)^6$$

bond-bond 与 bond-angle 交叉项与 class2 完全一致,无需重新定义。

4.2 系数定义

E_a需要设置 6 个系数:

系数单位说明
θ₀度平衡角度,内部转弧度
K₂能量energy/radian²
K₃能量energy/radian³
K₄能量energy/radian⁴
K₅能量energy/radian⁵
K₆能量energy/radian⁶

bb(M、r₁、r₂)与ba(N₁、N₂、r₁、r₂)系数的数量、顺序与 class2 完全相同,data 文件的BondBond Coeffs/BondAngle Coeffs节格式也一致。源码 src/MOFFF/angle_class2_p6.cpp 中coeff()默认分支校验narg != 7(含类型号 + 6 个系数),compute()中de_angle累加至6.0*k6*dtheta5,与六阶展开对应;extract()额外支持k5/k6关键字。该风格的单元测试见 unittest/force-styles/tests/angle-class2_p6.yaml。


五、系数配置的通用规则与注意事项

  1. 三种系数都必须设置:源码中coeff()在每轮设置后检查setflag_a、setflag_bb、setflag_ba三个标志,三者全置位才把setflag置 1,即角度主项、BondBond、BondAngle 缺一不可,否则该角度类型视为未完整定义。
  2. 通配符与覆盖:angle_coeff * ...对所有类型生效,之后可对单类型覆盖(angle_coeff 2 ...)。此机制与通用 angle_coeff 规则一致。
  3. data 文件格式差异:输入脚本中bb/ba分支必须带关键字,而 data 文件的BondBond Coeffs/BondAngle Coeffs节必须省略关键字,仅按“类型号 + 数值”逐行列出(可用 read_data 读取)。
  4. 角度顺序一致性:data 文件 Angles 节 I、J、K 顺序必须与r₁(i–j)、r₂(j–k)对应(见 2.5 节)。
  5. 重启文件:三类系数均写入 restart 文件,可经 read_restart 恢复(源码中write_restart/read_restart逐一序列化 θ₀、K₂、K₃、K₄ 与全部交叉项系数)。
  6. 单位体系:θ₀一律以度输入(内部转弧度),键长相关系数(r₁、r₂)与能量相关系数的实际量纲取决于当前 units 设置(如real、metal、lj等);class2xe的α为 1/距离量纲,需与单位制一致。

六、三个风格选型建议

  • 固定键长体系的 COMPASS/PCFF 建模:首选class2,配合angle_style class2所属 CLASS2 包中的bond_style class2、dihedral_style class2、improper_style class2、pair_style lj/class2/coul/long等风格可构建完整的 Class II 力场。
  • 需要模拟键拉伸解离、材料断裂(至断裂点):使用class2xe,键势换为bond_style morse,交叉项同步指数化;断裂后需叠加fix bond/react或fix bond/break。
  • 需要更高阶非谐角度响应:使用class2/p6(MOFFF 包),交叉项保持二次型即可。

所有风格均无默认系数(Default: none),必须通过angle_coeff或 data 文件显式给出三类系数后才能开始模拟。


七、相关命令与文档

  • 系数设置:angle_coeff(见 doc/src/angle_coeff.rst)
  • 数据文件:read_data(doc/src/read_data.rst)、read_restart(doc/src/read_restart.rst)
  • 力场背景:COMPASS(Sun 1998)、ClassII-xe(Kemppainen 2025),完整动机见 doc/src/Howto_bioFF.rst
  • 包构建:CLASS2、MOFFF 的安装参见 doc/src/Build_package.rst
  • 断裂后处理:fix bond/react(doc/src/fix_bond_react.rst)、fix bond/break(doc/src/fix_bond_break.rst)
  • 科研
  • 科学计算
  • 高性能计算

【免费下载链接】lammps

Public development project of the LAMMPS MD software package

项目地址:https://gitcode.com/gh_mirrors/la/lammps
点击查看免费下载

相关推荐

创作声明:本文部分内容由AI辅助生成(AIGC),仅供参考

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询