资讯动态

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

发布时间:2026/10/15 2:04:21 来源:尧图企业网站定制
科研科学计算高性能计算【免费下载链接】lammpsPublic 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 differentr_c通常应大于r_0即截断距离大于平衡键长这样才能保证势阱形态完整若r_c r_0键在未达到平衡距离前就被截断物理上虽然语法允许但不符合典型使用场景U_min应为正值其绝对值等于势阱深度平衡点能量-U_min。4. 与 harmonic、harmonic/shift 的差异对比LAMMPS 提供了三个互为姊妹关系的谐振键样式它们的势能形式对比如下键样式势能公式截断行为能量零点harmonicE K(r - r_0)²无硬截断可配合bond_style之外的全局截断平衡点E 0harmonic/shiftE U_min/(r_0-r_c)² · [(r-r_0)² - (r_c-r_0)²]无硬截断势能随距离持续变化平衡点E -U_minr_c处E 0harmonic/shift/cut与harmonic/shift相同r r_c时势能与力均强制归零平衡点E -U_minr_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换算成内部弹簧常数kdouble 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、r1r1即文档中的r_c三个数组U_min不再单独存储。通过 extract() 接口外部代码如单元测试或其他命令可以按k、r0、r1三个键名取出这些内部成员。5.2 力与能量计算主循环核心计算在 compute() 实现 中逐条遍历邻居列表中的键neighbor-bondlistrsq delx*delx dely*dely delz*delz; r sqrt(rsq); if (rr1[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/cutLAMMPS 会提示未知键样式。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.fourmolunittest/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 结构补充自己的回归用例。8. 小结bond_style harmonic/shift/cut通过一个常数平移项和硬截断把标准谐振键改造为平衡点能量为-U_min、截断点能量为零、超出截断距离力与能皆为零的键势。其使用要点可归纳为三个系数按U_min, r0, r_c顺序通过 bond_coeff 配置r0与r_c不能相等在r ≤ r_c区间内与K U_min/(r0-r_c)²的标准harmonic键力学等价需编译 EXTRA-MOLECULE 包可使用/omp、/kk加速变体底层实现bond_harmonic_shift_cut.cpp直接以截断判断 弹簧常数换算的方式落地势能公式并有 YAML 回归测试 提供数值级验证。对于需要在键断裂距离之外保证零相互作用、同时保留谐振弹性行为的模拟场景如聚合物拉伸、键断裂建模harmonic/shift/cut是一个简单可靠的选择。赞分享科研科学计算高性能计算【免费下载链接】lammpsPublic development project of the LAMMPS MD software package项目地址https://gitcode.com/gh_mirrors/la/lammps点击查看免费下载相关推荐LAMMPS bond_style harmonic/shift 移位谐振键势详解公式、系数、源码与加速变体LAMMPS bond_style harmonic/shift 移位谐振键势详解公式、系数、源码与加速变体 本篇技术指南围绕 LAMMPS 官方文档 bon科研科学计算高性能计算LAMMPS bond_style harmonic 详解谐振键势的参数配置、源码原理与加速变体LAMMPS bond_style harmonic 详解谐振键势的参数配置、源码原理与加速变体 bond_style harmonic 是 LAMMPS 中科研科学计算高性能计算SeaTunnel Hive 源连接器实战指南从单表读取到整库同步、Kerberos 认证与对象存储适配SeaTunnel Hive 源连接器实战指南从单表读取到整库同步、Kerberos 认证与对象存储适配 SeaTunnel 的 Hive 源连接器 con科研科学计算高性能计算上一篇Adobe Downloader完整指南3分钟掌握macOS上Adobe全家桶免费下载下一篇AzerothCore完整指南从零快速搭起自己的WotLK游戏服务器4步安装教程创作声明:本文部分内容由AI辅助生成(AIGC),仅供参考

读完文章,也想定制专属网站?

尧图设计师 24 小时内与您沟通定制方案

免费获取报价 →
↑