LAMMPS angle_style fourier/simple 傅里叶角度势详解：公式、参数、力计算与单元测试实证

• 科研
• 科学计算
• 高性能计算

angle_style fourier/simple是 LAMMPS 中一个基于单次余弦调制的角度势（E = K[1 + c·cos(nθ)]），它由 EXTRA-MOLECULE 包提供，支持kk、omp等加速变体。本文围绕官方文档 angle_fourier_simple.rst 展开，结合源码 angle_fourier_simple.cpp、angle_fourier_simple_kokkos.cpp 与单元测试 angle-fourier_simple.yaml，完整覆盖该角度风格的公式、系数含义、平衡角解析解、θ→0/π 奇点处理、KOKKOS/OMP 加速实现，并给出可直接运行的示例输入脚本与可复现的测试基准值。

一、命令语法、包依赖与使用前提

该命令的官方文档 angle_fourier_simple.rst 给出如下最简形式与示例：

angle_style fourier/simple angle_coeff 100.0 -1.0 1.0

并声明：

• 加速变体：fourier/simple/kk（KOKKOS 包）、fourier/simple/omp（OPENMP 包）；
• 包限制：仅当 LAMMPS 编译时启用了 EXTRA-MOLECULE 包才可用，构建方法见 Build package 文档；
• 默认值：none（即angle_style默认不选任何角度风格）；
• 相关命令：angle_coeff。

一个重要的形式勘误：文档示例angle_coeff 100.0 -1.0 1.0省略了角度类型号，但源码coeff()要求恰好 4 个参数（类型 + K + C + N），见 angle_fourier_simple.cpp：

void AngleFourierSimple::coeff(int narg, char **arg) { if (narg != 4) error->all(FLERR, "Incorrect args for angle coefficients" + utils::errorurl(21)); ... int ilo, ihi; utils::bounds(FLERR, arg[0], 1, atom->nangletypes, ilo, ihi, error); double k_one = utils::numeric(FLERR, arg[1], false, lmp); double C_one = utils::numeric(FLERR, arg[2], false, lmp); double N_one = utils::numeric(FLERR, arg[3], false, lmp);

因此实际可运行的形式必须带类型号（utils::expand_type支持通配范围如1*3与labelmap定义的标签，见 input.cpp）：

angle_style fourier/simple angle_coeff 1 100.0 -1.0 1.0

单元测试 in.fourmol 驱动的输入正是这种形式：

angle_coeff 1 75.0 -1.0 1.0 angle_coeff 2 45.0 1.0 2.0 angle_coeff 3 50.0 1.0 3.0 angle_coeff 4 100.0 -0.5 1.5

构建要求（前提说明，非仓库修改）：EXTRA-MOLECULE 依赖 MOLECULE 包（Packages_details.rst），典型配置为：

# CMake 构建 cmake -B build -DLAMMPS_PACKAGE_EXTRA-MOLECULE=on [其他选项...] # Make 构建 make yes-EXTRA-MOLECULE

二、势函数与三个系数

文档给出的角度势为：

E = K [ 1.0 + c cos( n θ ) ]

每个角度类型需通过angle_coeff指定三个系数（也可写在 data 文件的 "Angle Coeffs" 段，或 restart 文件中）：

源码中三个系数以double *k, *C, *N数组按角度类型索引存放（angle_fourier_simple.h），分配在allocate()中（L178-L189）。此外该类实现了extract()接口，支持通过compute ... extract fourier/simple k、... C、... N提取这三个内部数组（L327-L334）。

构造函数中还置位born_matrix_enable = 1并实现born_matrix()（L40-L46、L296-L321），意味着该风格可支持需要角度 Born 矩阵的功能（如声子类分析）；从源码结构看，这是少数实现该接口的角度风格之一。

三、能量与力的推导及奇点处理

compute()（angle_fourier_simple.cpp）对neighbor->anglelist中每个角度 (i1, i2, i3, type)：

1. 角度余弦：取两条键矢量 del1 = x[i1]−x[i2]、del2 = x[i3]−x[i2]，得

   c = del1·del2 / (r1 r2)，并裁剪到 [-1, 1] θ = acos(c)

2. 能量：eangle = k[type] * (1.0 + C[type] * cos(n θ))，与文档公式一致；eflag时累加，evflag时经基类ev_tally分摊能量与 6 项 virial（非newton on bond时按 1/3 分摊给三个原子）。

3. 力：对 θ 求导得 dE/dθ = −K·C·N·sin(nθ)，再经几何分解为标量

   a = K C N sin(nθ) / sin(θ)

   力作用方式为（源码注释明确说明"力是 θ 的函数，不依赖 r1、r2，只依赖角矢量"）：

   f1 = a·( c/r1² · del1 − 1/(r1 r2) · del2 ) # 作用于 i1 f3 = a·( c/r2² · del2 − 1/(r1 r2) · del1 ) # 作用于 i3 f2 = −(f1 + f3) # 作用于 i2

   从源码结构看，这保证三原子合力为零（f1+f2+f3=0）且对 i2 的净力矩为零（del1×f1 + del2×f3 = 0），即纯"扭转"约束力。

4. θ→0/π 奇点：a 中的sin(nθ)/sin(θ)在 θ→0 或 π 处分母趋于 0，代码用泰勒展开的极限形式替换（L123-L137）：

   // handle sin(n th)/sin(th) singulatiries if (fabs(c) - 1.0 > SMALL) { a = k[type] * C[type] * N[type] * sin(nth) / sin(th); } else { if (c >= 0.0) { term = 1.0 - c; sgn = 1.0; } else { term = 1.0 + c; sgn = (fmod(N[type], 2.0) == 0.0) ? -1.0 : 1.0; } a = N[type] + N[type] * (1.0 - N[type] * N[type]) * term / 3.0; a = k[type] * C[type] * N[type] * sgn * a; }

   其依据是极限 lim(θ→0) sin(nθ)/sin(θ) = n，以及 θ→π 时奇偶性：n 为偶数时该比值趋于 −n，为奇数时趋于 +n（由fmod(N, 2.0)判断）。SMALL = 0.0001为触发阈值。KOKKOS 变体逐字复现了同一分支逻辑（angle_fourier_simple_kokkos.cpp#L183-L197），保证 CPU/GPU 数值一致（除舍入外）。

5. 单点能量：single()（L268-L292）对给定角度按最小镜像约定重算能量，供蒙特卡洛类修复等外部场景调用，公式与compute()完全相同。

四、平衡角度：θ_eq = π/n

类实现了equilibrium_angle()：

double AngleFourierSimple::equilibrium_angle(int i) { return (MY_PI / N[i]); }

（angle_fourier_simple.cpp#L221-L224）。当 c = −1 时能量在 cos(nθ) = +1（θ = π/n）处取最小值 K(1−1) = 0；当 c = +1 时 π/n 为极大、0 与 π 为极小。该解析解被单元测试直接校验（angle-fourier_simple.yaml）：

四类全部等于 π/n，印证了该实现与公式定义的一致性。

五、加速变体（kk / omp）的实现细节

文档声明fourier/simple/kk、fourier/simple/omp两种加速风格，并引用 accel_styles.rst 的通用说明：带后缀风格与基础风格功能相同，仅当对应包（KOKKOS、OPENMP 等）被编译启用时可用，可通过显式写后缀、-suffix命令行开关或suffix输入命令启用。

• KOKKOS 变体：类AngleFourierSimpleKokkos<DeviceType>继承自AngleFourierSimple（angle_fourier_simple_kokkos.h#L16-L18），kokkosable = 1；系数 k/C/N 以 dual-view 同步到设备端（allocate()中k_k/k_C/k_N，kokkos.cpp#L244-L257）；coeff()与read_restart()在写入主机后显式modify_host()（#L263-L280、#L286-L301）；
• compute()用Kokkos::parallel_reduce / parallel_for在 0..nanglelist 上并行，按newton_bond与evflag分 4 个模板实例（#L103-L117），力以Kokkos::Atomic原子累加写入f，避免非牛顿模式下的竞争写；
• 能量与 virial 的原子分摊逻辑与 CPU 版一致（ev_tally，#L312-L382）；
• 注释说明 MC 类修复会在run_style verlet/kk之外调用该风格，因此compute()内自行完成数据同步与modified标记（#L87-L91）；
• OMP 变体：类AngleFourierSimpleOMP继承AngleFourierSimple与ThrOMP，仅重写compute()（angle_fourier_simple_omp.h#L32-L41），系数处理复用基类。

从源码结构看，两个加速变体都没有重写single()、born_matrix()等接口，相关功能仍走基类 CPU 路径。

六、完整可运行示例

以下脚本基于单元测试使用的数据文件 data.fourmol（含 5 原子类型、30 个角度、4 种角度类型），结构参照 in.fourmol：

variable newton_pair index on variable newton_bond index on variable bond_factor index 0.10 variable angle_factor index 0.25 variable dihedral_factor index 0.50 atom_style full atom_modify map array neigh_modify delay 2 every 2 check no units real timestep 0.1 newton ${newton_pair} ${newton_bond} special_bonds lj/coul ${bond_factor} ${angle_factor} ${dihedral_factor} pair_style zero 8.0 bond_style zero angle_style fourier/simple # 或 fourier/simple/kk（需 -suffix / suffix 命令） dihedral_style zero improper_style zero read_data data.fourmol # 需含 Angle Coeffs 段，或读后用 angle_coeff 设置 angle_coeff 1 75.0 -1.0 1.0 angle_coeff 2 45.0 1.0 2.0 angle_coeff 3 50.0 1.0 3.0 angle_coeff 4 100.0 -0.5 1.5 compute 1 all energy/angle compute 2 all stress/atom dump 1 all custom 100 traj.lammps id type x y z fx fy fz run 100

使用要点：

• 若 data 文件不含 "Angle Coeffs" 段，必须在angle_style之后用angle_coeff逐一设置；系数行格式为type K c n（angle_coeff.rst 说明 data 文件写法与命令参数相同，只是不允许通配符）；
• angle_style必须在任何angle_coeff之前出现，且需先有模拟盒子（read_data/create_box），见 input.cpp#L1370-L1383 的检查顺序；
• atom_style必须允许角度（测试用full，prerequisites 中声明为atom full+angle fourier/simple）；
• 使用kk变体时，系数在coeff()中自动同步到 GPU 设备视图，无需额外命令。

七、单元测试：可复现的数值基准

单元测试 angle-fourier_simple.yaml 为该风格提供了回归基准（epsilon: 2.5e-13，即 1e-13 量级对拍容差），关键断言值：

• natoms: 29，29 个原子的逐原子初始力init_forces（全列表见 yaml 文件）；
• init_energy: 2474.0748013590646，6 项初始应力init_stress；
• 运行后run_energy: 2473.1257456970907及run_stress、run_forces；
• equilibrium四个平衡角（见上表）。

这些数值可直接用于验证自编译二进制（含kk/omp变体）的数值一致性：在相同输入下重放该测试，力与能量应落在 epsilon 之内。

八、参考文件索引

• 科研
• 科学计算
• 高性能计算