LAMMPS bond_style harmonic/shift/cut 详解：位移截断谐振键势的参数、原理与实现-育师

科研
科学计算
高性能计算

【免费下载链接】lammps

Public development project of the LAMMPS MD software package

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

导读

bond_style harmonic/shift/cut是 LAMMPS 中一种将谐振键势整体向下平移、并在临界距离r_c处硬截断的键模型，其势能曲线在平衡距离r_0处取最小值-U_min，在r_c处恰好为零，超出r_c后势能与力同时归零。本文以官方文档 doc/src/bond_harmonic_shift_cut.rst 为主体，结合仓库源码与单元测试，完整讲解该键样式的势能公式、三个系数（U_min、r_0、r_c）的物理含义与配置方法、与harmonic及harmonic/shift的差异、加速变体（/omp、/kk）的使用方式，以及它从参数解析到力/能量计算的底层实现细节。读完本文，你将能在自己的输入脚本中正确配置并高效使用这一键势，并理解其与标准谐振势的换算关系。

1. 命令语法与快速上手

bond_style harmonic/shift/cut在输入脚本中的使用方式非常简单，语法如下：

bond_style harmonic/shift/cut

对应的官方示例：

bond_style harmonic/shift/cut bond_coeff 5 10.0 0.5 1.0

其中bond_coeff的行为第 5 号键类型依次指定三个系数（详见第 3 节）：U_min = 10.0（能量单位）、r_0 = 0.5（距离单位）、r_c = 1.0（距离单位）。

该键样式没有额外的默认参数（官方文档的 Default 一节为 none），所有行为完全由bond_coeff给定的三个系数决定。要使用它，只需保证 LAMMPS 在编译时启用了EXTRA-MOLECULE包（限制条件详见第 6 节）。

1.1 完整的可运行输入片段

将上面的命令放入一个最小的输入脚本中即可运行，例如：

units real atom_style full read_data data.fourmol bond_style harmonic/shift/cut bond_coeff 1 250.0 1.5 2.0 bond_coeff 2 300.0 1.1 1.5 bond_coeff 3 350.0 1.3 1.5 run 100

需要说明的是，bond_coeff既可以在输入脚本中通过命令逐类型设置，也可以在数据文件（由 read_data 读取）或重启文件（由 read_restart 读取）中携带对应格式的系数，这一点与 LAMMPS 其他键样式完全一致。

2. 势能公式与物理含义

harmonic/shift/cut键势的势能表达式为：

$$ E = \frac{U_{\text{min}}}{(r_0 - r_c)^2} \left[ (r - r_0)^2 - (r_c - r_0)^2 \right] $$

其中：

r_0：平衡键长（equilibrium bond distance），即势能取极值的位置；
r_c：临界距离（critical distance），即势能截断位置；
U_min：势阱深度，势能在r_0处的取值为-U_min。

该势能的几个关键性质：

平衡点势能为负：在r = r_0处，(r - r_0)^2 = 0，代入公式可得E = -U_min，即势能曲线被整体向下平移，最小值不再是 0；
截断点势能为零：在r = r_c处，(r - r_0)^2 = (r_c - r_0)^2，两式相消得E = 0；
超过截断距离完全失效：对于r > r_c的键对，势能为零，对应的力也为零。这一点正是harmonic/shift/cut与harmonic/shift（其势能曲线随距离无限延伸）的本质区别。

值得注意的是，harmonic/shift/cut与标准 bond_style harmonic 的势能表达式

$$ E = K (r - r_0)^2 $$

相比，多了一个常数项-U_min（因为U_min / (r_0 - r_c)^2 × (r_c - r_0)^2 = U_min）。由于力的计算只对势能求导，常数项不影响任何力，两者的力学行为在r ≤ r_c范围内完全等价，差异仅体现在势能的绝对数值（能量标签）上。

2.1 与标准谐振势的换算关系

官方文档明确指出，若想得到与harmonic/shift/cut在r ≤ r_c范围内相同的力学行为，标准harmonic键势的等效弹簧常数可按下式换算：

$$ K = \frac{U_{\text{min}}}{(r_0 - r_c)^2} $$

换言之，bond_coeff 5 10.0 0.5 1.0配置下的harmonic/shift/cut，与bond_coeff 5 40.0 0.5（即K = 10.0 / (0.5 - 1.0)² = 40.0）的标准harmonic键在r ≤ 1.0范围内给出相同的力。这一换算关系也在源码的参数处理中得到印证（见第 5 节）。

3. 三个系数的定义与取值要点

通过bond_coeff命令，每个键类型必须定义以下三个系数（顺序固定）：

系数	符号	物理意义	量纲
第 1 个	`U_min`	势能最小值深度（在`r_0`处势能为`-U_min`）	能量（如 kcal/mol、eV）
第 2 个	`r_0`	平衡键距离（势能极小值位置）	距离（如 Å）
第 3 个	`r_c`	临界/截断距离（势能为零并开始截断的位置）	距离（如 Å）

取值要点（可由源码的校验逻辑推导，见第 5.1 节）：

r_0与r_c不能相等。若两者相同，(r_0 - r_c)² = 0，等效弹簧常数K将趋于无穷大，公式失去意义，LAMMPS 会直接报错终止（源码中对应错误信息为Bond harmonic/shift/cut r0 and r1 must be different）；
r_c通常应大于r_0，即截断距离大于平衡键长，这样才能保证势阱形态完整；若r_c < r_0，键在未达到平衡距离前就被截断，物理上虽然语法允许，但不符合典型使用场景；
U_min应为正值，其绝对值等于势阱深度（平衡点能量-U_min）。

4. 与 harmonic、harmonic/shift 的差异对比

LAMMPS 提供了三个互为姊妹关系的谐振键样式，它们的势能形式对比如下：

键样式	势能公式	截断行为	能量零点
`harmonic`	`E = K(r - r_0)²`	无硬截断（可配合`bond_style`之外的全局截断）	平衡点`E = 0`
`harmonic/shift`	`E = U_min/(r_0-r_c)² · [(r-r_0)² - (r_c-r_0)²]`	无硬截断，势能随距离持续变化	平衡点`E = -U_min`，`r_c`处`E = 0`
`harmonic/shift/cut`	与`harmonic/shift`相同	`r > r_c`时势能与力均强制归零	平衡点`E = -U_min`，`r_c`处`E = 0`

三者中，harmonic依赖 MOLECULE 包，而harmonic/shift与harmonic/shift/cut同属 EXTRA-MOLECULE 包。官方文档的 Related commands 部分同时列出了 bond_coeff、delete_bonds、bond_harmonic 与 bond_style harmonic/shift 作为关联命令，方便读者对照查阅。

harmonic/shift/cut的典型适用场景是：希望在截断距离之外严格不产生任何相互作用（力与能量都为零），从而避免长键在拉伸断裂后仍残留非零力导致数值问题，同时又要保持r ≤ r_c区间内标准的谐振力学行为。

5. 源码级实现解析

harmonic/shift/cut的核心实现位于 EXTRA-MOLECULE 包中，由头文件 bond_harmonic_shift_cut.h 与实现文件 bond_harmonic_shift_cut.cpp 构成（贡献作者为 Carsten Svaneborg）。下面沿参数解析 → 力/能量计算 → 辅助接口三条线索展开。

5.1 参数解析：从 U_min、r0、r1 到等效弹簧常数

在 coeff() 实现中，bond_coeff要求恰好 4 个参数（类型范围 + 3 个系数），否则报错。核心逻辑是先把U_min换算成内部弹簧常数k：

double Umin = utils::numeric(FLERR,arg[1],false,lmp); // energy at minimum double r0_one = utils::numeric(FLERR,arg[2],false,lmp); // position of minimum double r1_one = utils::numeric(FLERR,arg[3],false,lmp); // position where energy = 0 = cutoff if (r0_one == r1_one) error->all(FLERR,"Bond harmonic/shift/cut r0 and r1 must be different"); for (int i = ilo; i <= ihi; i++) { k[i] = Umin/((r0_one-r1_one)*(r0_one-r1_one)); r0[i] = r0_one; r1[i] = r1_one; setflag[i] = 1; }

这正是文档中K = U_min / [(r_0 - r_c)²]公式的代码形态：内部只保存k、r0、r1（r1即文档中的r_c）三个数组，U_min不再单独存储。通过 extract() 接口，外部代码（如单元测试或其他命令）可以按"k"、"r0"、"r1"三个键名取出这些内部成员。

5.2 力与能量计算主循环

核心计算在 compute() 实现中，逐条遍历邻居列表中的键（neighbor->bondlist）：

rsq = delx*delx + dely*dely + delz*delz; r = sqrt(rsq); if (r>r1[type]) continue; // 硬截断：超出 r_c 后势能与力都为零 dr = r - r0[type]; rk = k[type] * dr; if (r > 0.0) fbond = -2.0*rk/r; else fbond = 0.0; if (eflag) ebond = k[type]*(dr*dr -(r0[type]-r1[type])*(r0[type]-r1[type]));

这段代码与势能公式一一对应：

截断判定：r > r1[type]（即r > r_c）时直接continue，跳过该键，不贡献任何力与能量——这就是文档所述"势能及其力在r > r_c时为零"的代码实现；
力表达式：F = -dE/dr = -2k(r - r0)，程序按单位矢量展开为fbond = -2·k·(r-r0)/r，再沿delx/dely/delz分量分别施加到两个原子（并遵循newton_bond与局部原子归属的判断）；
能量表达式：ebond = k·[(r-r0)² - (r0-r_c)²]，与文档公式完全一致，r0-r1即r0-r_c。

5.3 配套辅助接口

除主循环外，该类还实现了若干供 LAMMPS 框架与后处理使用的接口：

equilibrium_distance(int i)：返回r0[i]，即平衡键长，供依赖平衡几何的模块（如键长约束初始化）查询；
single(type, rsq, ...)：对单个键对独立计算力与能量（同样先做r > r1[type]的截断判断），供如compute bond/local等按键输出属性的命令使用；
born_matrix(type, rsq, ...)：提供 Born 矩阵相关的一阶、二阶导数，du2 = 2k，供弹性性质计算使用；
write_restart / read_restart：将k、r0、r1三个数组写入重启文件并在读取后广播，保证重启后系数一致；
write_data：将系数写回数据文件时，回写的是k·(r0-r1)²即原始U_min，保证数据文件中的系数格式与输入脚本一致（%d %g %g %g对应类型号、U_min、r0、r1）。

5.4 加速变体的实现

官方文档声明该样式带有harmonic/shift/cut/omp与harmonic/shift/cut/kk两个加速变体（Accelerator Variants），它们功能等价，仅在数值精度（舍入）上可能有微小差异：

OMP 变体：bond_harmonic_shift_cut_omp.cpp 通过loop_setup_thr将键列表按线程切分，在#pragma omp parallel区域内用模板参数EVFLAG/EFLAG/NEWTON_BOND组合展开同一套力/能量计算，再经reduce_thr归并全局累加量；
KOKKOS 变体：bond_harmonic_shift_cut_kokkos.cpp 将k/r0/r1以双视图（host/device）形式管理，键循环通过Kokkos::parallel_for或parallel_reduce在设备端并行执行，并实现了ev_tally的能量/维里累加逻辑。

两个加速变体分别隶属于 OPENMP 与 KOKKOS 包，只有在编译时启用对应包才可用（详见 Build package 文档）。启用方式有三种：在输入脚本中显式书写带后缀的样式名、使用 LAMMPS 命令行-suffix开关、或在脚本中用suffix命令统一追加后缀（文档 accel_styles.rst 对此有完整说明）。

6. 使用限制：EXTRA-MOLECULE 包

官方文档的 Restrictions 一节明确指出：harmonic/shift/cut只能在编译时启用了 EXTRA-MOLECULE 包的 LAMMPS 构建中使用。启用方式参考 Build package 文档：

传统 make 构建：在src目录下执行make yes-extra-molecule，然后重新编译（make构建时对包目录执行make yes-name/make no-name即可安装/卸载对应包，安装后必须重新 build 才生效）；
CMake 构建：在配置阶段启用 EXTRA-MOLECULE 选项（或直接使用包含该包的预设）。

可以在启动 LAMMPS 前用-h帮助输出确认当前构建是否包含该样式。若未启用该包而直接在输入脚本中使用bond_style harmonic/shift/cut，LAMMPS 会提示未知键样式。

7. 单元测试与验证

仓库的回归测试体系为harmonic/shift/cut提供了完整的数值验证，测试配置位于 unittest/force-styles/tests/bond-harmonic_shift_cut.yaml，由 test_bond_style.cpp 驱动（测试通过 CMake 的add_test(NAME ... COMMAND test_bond_style ...)机制注册，见 unittest/force-styles/CMakeLists.txt）。

该测试文件的关键信息：

lammps_version: 12 Jun 2025 epsilon: 2.5e-13 prerequisites: | atom full bond harmonic/shift/cut input_file: in.fourmol bond_style: harmonic/shift/cut bond_coeff: | 1 250.0 1.5 2.0 2 300.0 1.1 1.5 3 350.0 1.3 1.5 4 650.0 1.2 0.5 5 450.0 1.0 1.3 equilibrium: 5 1.5 1.1 1.3 1.2 1 extract: | k 1 r0 1 r1 1

从中可以提取几个有价值的实操细节：

多键类型批量配置：示例一口气为 5 种键类型配置了不同的U_min / r0 / r1，其中类型 4 的r_c = 0.5 < r_0 = 1.2，说明语法层面允许r_c < r_0（但物理上并不常见）；
平衡键长验证：equilibrium字段列出5种类型的平衡距离分别为1.5、1.1、1.3、1.2、1.0，正是各r0值——与源码中equilibrium_distance()返回r0[i]的行为一致；
内部成员导出：extract字段要求按k、r0、r1导出内部数组，验证了 extract() 接口的契约；
精度标准：测试容差epsilon = 2.5e-13，对能量、应力与逐原子力做数值比对，测试输入in.fourmol（unittest/force-styles/tests/in.fourmol）采用units real与atom_style full，并配合special_bonds lj/coul处理 1-2/1-3/1-4 排除规则。

此外，测试框架支持newton on/off与newton_bond on/off的组合验证（in.fourmol中通过变量newton_pair、newton_bond控制），确保harmonic/shift/cut在牛顿第三定律开关两种模式下都给出一致结果。如果你修改或自定义了该键样式，可参照此 YAML 结构补充自己的回归用例。