LAMMPS OPENMP 包风格开发指南：从 PairLJCutOMP 到 ThrOMP 多线程机制详解-平芜编程栈

科研
科学计算
高性能计算

【免费下载链接】lammps

Public development project of the LAMMPS MD software package

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

本文基于 LAMMPS 官方开发者文档doc/src/Developer_write_openmp.rst，结合 src/OPENMP 目录下的真实源码，完整讲解如何为 LAMMPS 的 OPENMP 包编写一个新的多线程加速风格（style）：从继承体系、compute()改写、ThrData线程本地缓存、模板化求力内核，到 EAM/SW/Tersoff 等多体势的特殊移植难点与 bonded 风格的移植套路。读完本文，你能够独立理解任意一个*_omp风格的源码结构，并掌握将新的 pair/bond/angle 风格移植到 OPENMP 包的全部关键步骤。

OPENMP 包定位与构建系统要求

OPENMP 包是 LAMMPS 的"加速器包"（accelerator package）之一，它通过 OpenMP pragma 以多线程方式加速大量 LAMMPS 风格，与基于 Kokkos 的加速器包遵循的并行化策略是正交的。高层原理可参阅开发者手册中的 OpenMP Parallelism 章节。

将新风格加入 OPENMP 包时有两条硬性约束：

代码位置：新的*_omp风格源文件必须放在src/OPENMP目录中。该目录当前包含数百个多线程风格，命名规律清晰：pair_lj_cut_omp.h/.cpp、bond_harmonic_omp.h/.cpp、angle_cosine_omp.h/.cpp、dihedral_charmm_omp.h/.cpp、improper_ring_omp.h/.cpp等，此外还有fix_nh_omp、fix_npt_omp、pppm_omp、npair_nsq_omp等热力学、长程和邻居列表的加速实现。
双构建系统支持：贡献代码需要同时支持 CMake 构建流程和传统的 GNU make 流程。在 CMake 侧，cmake/CMakeLists.txt中有set(SUFFIX_PACKAGES CORESHELL GPU KOKKOS OPT INTEL OPENMP)，OPENMP 包通过LAMMPS_PKG_OMP选项启用并自动开启 OpenMP 编译支持；在 make 流程下则使用src/MAKE/Makefile.omp（其中定义了-DOPENMP -fopenmp等开关）。

案例研究：成对加和势 pair_lj_cut

对许多 pair 风格而言，添加 OpenMP 支持意味着创建一个继承自串行风格的派生类。下面以pair_lj_cut及其加速版pair_lj_cut/omp为例，逐段对照源码讲解。

头文件：双重继承 PairLJCut 与 ThrOMP

openmp 版头文件的骨架与文档描述一致（略有扩充，文档示例做了精简）：

#ifndef LMP_PAIR_LJ_CUT_OMP_H #define LMP_PAIR_LJ_CUT_OMP_H #include "pair_lj_cut.h" #include "thr_omp.h" namespace LAMMPS_NS { class PairLJCutOMP : public PairLJCut, public ThrOMP { public: PairLJCutOMP(class LAMMPS *); void compute(int, int) override; double memory_usage() override; private: template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void eval(int ifrom, int ito, ThrData *const thr); }; } // namespace LAMMPS_NS #endif

要点：

类继承自PairLJCut（串行基类）和ThrOMP（线程工具基类）。ThrOMP负责管理 OpenMP 专属数据，例如线程本地的力与能量累加数组，定义于 thr_omp.h。
只 override 了两个方法：compute()（真正的求力入口）和memory_usage()（报告线程本地缓冲的额外内存）。其余方法全部继承。
私有模板函数eval<EVFLAG, EFLAG, NEWTON_PAIR>是实际的热核循环。

风格注册发生在头文件顶部的#ifdef PAIR_CLASS分支中：PairStyle(lj/cut/omp,PairLJCutOMP);。LAMMPS 的pair_style lj/cut/omp命令因此能找到PairLJCutOMP类。

构造函数与"为什么只重写少数方法"

真实的 pair_lj_cut_omp.cpp 构造函数：

PairLJCutOMP::PairLJCutOMP(LAMMPS *lmp) : PairLJCut(lmp), ThrOMP(lmp, THR_PAIR) { suffix_flag |= Suffix::OMP; respa_enable = 0; cut_respa = nullptr; }

THR_PAIR参数告诉ThrOMP基类这是 pair 风格的加速。查看 thr_omp.h 中的枚举：THR_NONE、THR_PAIR、THR_BOND、THR_ANGLE、THR_DIHEDRAL、THR_IMPROPER、THR_KSPACE、THR_CHARMM、THR_FIX、THR_INTGR，它直接影响ev_setup_thr()为哪一类风格初始化线程本地数组。
respa_enable = 0是 OMP 加速风格的常见做法——RESPA 加速通常需要单独的实现逻辑（cut_respa = nullptr则显式禁用 RESPA 切换路径）。
suffix_flag |= Suffix::OMP标记该风格使用/omp后缀，这是 LAMMPS 的 suffix 机制（声明于 suffix.h），保证带后缀的风格名在注册表中与串行风格区分开。

大多数方法（coeff()、init_one()、write_restart()等）直接从PairLJCut继承，原因有三：

串行类已经处理了全部数据结构和输入解析；
消耗绝大部分 CPU 时间、且每步都会被调用的方法是compute()，只需加速它；
single()没有被 override，因为它无法被 OpenMP 高效并行化——它只是compute()主循环最内层的一次迭代。

线程工具函数：loop_setup_thr、get_thr、reduce_thr

ThrOMP基类和 LAMMPS 基础设施提供了几个管理并行执行的关键方法：

loop_setup_thr(int inum, int nthreads)计算每个线程处理的原子（或邻居）范围，实现就在 thr_omp.h 中：

inline void loop_setup_thr(int &ifrom, int &ito, int &tid, int inum, int OMP_ONLY_PARAM(nthreads)) { #if defined(_OPENMP) tid = omp_get_thread_num(); // each thread works on a fixed chunk of atoms. const int idelta = 1 + inum / nthreads; ifrom = tid * idelta; ito = ((ifrom + idelta) > inum) ? inum : ifrom + idelta; #else tid = 0; ifrom = 0; ito = inum; #endif }

即"静态分块"策略：每个线程负责固定的连续区间，非 OpenMP 编译时退化为单线程全区间，保证了同一份代码两种构建方式都能工作。

FixOMP::get_thr(int tid)按线程号取ThrData对象，见 fix_omp.h：

ThrData *get_thr(int tid) { return thr[tid]; }

ThrData（定义于 thr_data.h）保存每个线程私有的力、能量、virial 缓冲区，使多线程在并行区内各自累加贡献而无锁竞争。

reduce_thr(ThrData* thr, int eflag, int vflag)在并行区结束后把各线程的力、能量、virial 归约到全局数组。thr_omp.cpp 中的实现细节值得注意：

全局量（总能量、pair virial 等）通过#pragma omp critical串行归约，开销小；
逐原子量（eatom、vatom等）通过data_reduce_thr()在多个线程间并行归约，方式与力归约相同——这正是"无锁并行累加 + 事后归约"策略的体现。

ev_setup_thr 与 ev_tally_thr：线程本地化的能量/virial 统计

这两个方法替代了非 OpenMP 的对应函数ev_setup()和ev_tally()：

ev_setup_thr(eflag, vflag, nall, eatom, vatom, cvatom, thr)：根据thr_style的位标志，把全局 per-atom 数组按eatom + tid*nall偏移映射到当前线程的切片并清零，见 thr_omp.cpp 中case THR_PAIR分支：

if (thr_style & THR_PAIR) { if (eflag & ENERGY_ATOM) { thr->eatom_pair = eatom + tid*nall; if (nall > 0) memset(&(thr->eatom_pair[0]),0,nall*sizeof(double)); } ... }

也就是说 per-atom 能量/virial 的底层大数组被切成nthreads份，每线程独占一份，天然无竞争。

ev_tally_thr(i, j, nlocal, newton, evdwl, ecoul, fpair, delx, dely, delz, thr)：不直接更新全局数组，而是更新ThrData内的线程本地累加器。Pair 版实现中先调用e_tally_thr()（能量）与v_tally_thr()（virial，用delx*delx*fpair等 6 个分量累加），最后若注册了pair_tally回调，则用#pragma omp critical保护回调调用——因为回调本身不是线程安全的。

此外 thr_omp.h 还提供了check_error_thr()（线程安全的报错机制：任意线程置位thr_error，由线程 0 统一调用error->one()终止）以及dbl3_t/dbl4_t等结构体别名，后者用于把double**的atom->x[0]解释为结构体指针，帮助编译器做更好的向量化优化。

compute() 与模板化 eval() 内核

PairLJCutOMP::compute() 完整流程：

void PairLJCutOMP::compute(int eflag, int vflag) { ev_init(eflag,vflag); const int nall = atom->nlocal + atom->nghost; const int nthreads = comm->nthreads; const int inum = list->inum; #if defined(_OPENMP) #pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag,vflag) #endif { int ifrom, ito, tid; loop_setup_thr(ifrom, ito, tid, inum, nthreads); ThrData *thr = fix->get_thr(tid); thr->timer(Timer::START); ev_setup_thr(eflag, vflag, nall, eatom, vatom, nullptr, thr); if (evflag) { if (eflag) { if (force->newton_pair) eval<1,1,1>(ifrom, ito, thr); else eval<1,1,0>(ifrom, ito, thr); } else { if (force->newton_pair) eval<1,0,1>(ifrom, ito, thr); else eval<1,0,0>(ifrom, ito, thr); } } else { if (force->newton_pair) eval<0,0,1>(ifrom, ito, thr); else eval<0,0,0>(ifrom, ito, thr); } thr->timer(Timer::PAIR); reduce_thr(this, eflag, vflag, thr); } // end of omp parallel region }

说明（注意源码中的小细节）：

并行区用#pragma omp parallel开启；LMP_DEFAULT_NONE/LMP_SHARED(...)是omp_compat.h里对 OpenMP 版本差异做兼容的宏，最终效果等价于default(none) shared(eflag,vflag)。
ev_init(eflag,vflag)根据命令行的energy/compute pressure等设置置位evflag，之后每个线程据此选择模板实例。
eval<EVFLAG,EFLAG,NEWTON_PAIR>用模板参数替代了热循环里的运行时条件判断。上例 8 种模板实例覆盖了 evflag×eflag×newton_pair 的全部组合，分支在编译期消除——这正是文档所说的"Branch Elimination"优化。
注意并行区结尾：reduce_thr()也在#pragma omp parallel区块内部调用，各线程先归约自己那份数据，再在函数内部的 barrier/critical 上同步。

再看 eval() 热核的关键片段：

const auto * _noalias const x = (dbl3_t *) atom->x[0]; auto * _noalias const f = (dbl3_t *) thr->get_f()[0]; const int * _noalias const type = atom->type; const int * _noalias const ilist = list->ilist; ... for (int ii = iifrom; ii < iito; ++ii) { const int i = ilist[ii]; ... for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; factor_lj = special_lj[sbmask(j)]; j &= NEIGHMASK; ... if (rsq < cutsqi[jtype]) { r2inv = 1.0/rsq; r6inv = r2inv*r2inv*r2inv; forcelj = r6inv * (lj1i[jtype]*r6inv - lj2i[jtype]); fpair = factor_lj*forcelj*r2inv; fxtmp += delx*fpair; ... if (NEWTON_PAIR || j < nlocal) { f[j].x -= delx*fpair; ... } if (EFLAG) { evdwl = r6inv*(lj3i[jtype]*r6inv-lj4i[jtype]) - offseti[jtype]; evdwl *= factor_lj; } if (EVFLAG) ev_tally_thr(this,i,j,nlocal,NEWTON_PAIR, evdwl,0.0,fpair,delx,dely,delz,thr); } } f[i].x += fxtmp; ... }

两个性能要点与文档完全对应：

模板消除分支：if (NEWTON_PAIR || j < nlocal)、if (EFLAG)等条件中的标志是模板参数，编译器为每个组合生成专门的机器码，热循环里不再有可翻转的分支，对分支预测器非常友好。
_noalias提示：告诉编译器各指针（如x、f）不指向重叠内存，允许激进地重排 load/store 和展开循环。
原子i的力先在寄存器级临时变量fxtmp/fytmp/fztmp中累加、循环结束后一次性写回f[i]；原子j的力直接累加进线程本地f数组（thr->get_f()返回的是本线程切片）。
special_lj[sbmask(j)]保留了对bond/angle/coefficient special等特殊对系数的支持，与串行版行为一致。

memory_usage()则把线程本地缓冲的内存也计入统计：memory_usage_thr() + PairLJCut::memory_usage()，保证print等命令报告的值准确。

移植复杂多体势：PairEAMOMP、PairSWOMP、PairTersoffOMP

多体势（Stillinger-Weber、Tersoff 等）与成对加和势（LJ）有本质区别，移植时需注意三类差异：

邻居列表要求：成对加和势通常用"半"邻居列表即可；多体势一般需要"full"邻居列表，因为计算三体项i-j-k时既要遍历i的邻居，也要遍历j的邻居。
计算结构：不是简单的嵌套双层循环，而是涉及"邻居的邻居"的复杂内核，并行分区与数据局部性要求更高。
通信：由于三体项需要访问j截断半径内所有原子的坐标（即使j本身是 ghost 原子），需要同步的 ghost 数据量更大。

以 EAM 为例，文档提到多遍扫描（密度累加一遍、求力一遍）。当前源码 pair_eam_omp.cpp 的实际组织方式值得细读——它并不是把两个并行区简单串联，而是在一个#pragma omp parallel区内用多个串行 pass 加sync_threads()（即#pragma omp barrier）编排，且 MPI 通信只由#pragma omp master主线程执行：

// pass 1: 各线程并行累加本线程切片上的原子密度 for (ii = iifrom; ii < iito; ii++) { ... rho_t[i] += ...; ... } sync_threads(); // 所有线程完成密度累加 if (NEWTON_PAIR) { data_reduce_thr(rho, nall, nthreads, 1, tid); // 多线程并行归约密度 sync_threads(); #pragma omp master { comm->reverse_comm(this); } // 仅主线程做 MPI 反向通信 sync_threads(); } // pass 2: 计算嵌入能导数 fp、能量 phi（含超线性外推警告） for (ii = iifrom; ii < iito; ii++) { embedding_index<HE>(rho[i],m,p); ... } sync_threads(); #pragma omp master { comm->forward_comm(this); } // 把 fp 同步到 ghost 原子 sync_threads(); // pass 3: 各线程并行求力 for (ii = iifrom; ii < iito; ii++) { ... fpair = -scale_i[jtype]*psip*recip; ... }

模板参数扩展为eval<EVFLAG, EFLAG, NEWTON_PAIR, HE>，其中HE表示eam/he变体的表边界处理。文档中"两遍并行 + 中间 forward_comm"的示意代码是对这一编排的抽象简化；从当前源码看，实际的同步点是sync_threads()屏障、data_reduce_thr()并行归约，以及由主线程执行的comm->reverse_comm(this)/comm->forward_comm(this)。

对PairSWOMP和PairTersoffOMP（pair_sw_omp.cpp、pair_tersoff_omp.cpp），加速思路是并行化嵌套循环同时保持线程安全：

数据分区：compute()必须把（full）邻居列表在线程间仔细分区，保证每个线程处理互不重叠的相互作用集；
Tallying：与成对风格同样用ev_tally_thr()，但多体势中 2 体项与 3 体项的贡献常需分别累积——ThrOMP为此提供了专门的线程安全统计函数，如ev_tally3_thr()（注释标明"called by SW and hbond potentials"，按drji*fj + drki*fk计算三体 virial，能量按三分之一分摊）和ev_tally4_thr()（AIREBO 的四体项，能量按四分之一分摊）；
内核特化：eval()依旧重度模板化，处理 Newton 第三定律开/关、能量开关等场景。

移植的核心挑战在于确保多体逻辑在多线程环境下的正确性：对共享邻居列表和线程本地原子数据，要通过ThrData结构严格划分所有权，屏障只出现在真正的全局依赖点（如 EAM 的密度归约、ghost 同步）。

加速 bonded 风格：bond/angle/dihedral/improper 的通用套路

bonded 风格与 pair 风格遵循完全相同的哲学：BondHarmonicOMP继承BondHarmonic与ThrOMP（用THR_BOND标志），同样只 overridecompute()。以 bond_harmonic_omp.cpp 为例，结构与 pair 版几乎一一对应：

void BondHarmonicOMP::compute(int eflag, int vflag) { ev_init(eflag,vflag); const int nall = atom->nlocal + atom->nghost; const int nthreads = comm->nthreads; const int inum = neighbor->nbondlist; // 注意：分区对象是键列表长度 #if defined(_OPENMP) #pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag,vflag) #endif { int ifrom, ito, tid; loop_setup_thr(ifrom, ito, tid, inum, nthreads); ThrData *thr = fix->get_thr(tid); ... if (inum > 0) { if (evflag) { if (eflag) { if (force->newton_bond) eval<1,1,1>(ifrom, ito, thr); else eval<1,1,0>(ifrom, ito, thr); } else { ... } } else { ... } } thr->timer(Timer::BOND); reduce_thr(this, eflag, vflag, thr); } }

与 pair 风格仅有的差异：

分区的inum取自neighbor->nbondlist（键的数量），而非邻居列表条目数；
模板参数是<EVFLAG, EFLAG, NEWTON_BOND>；
eval()中键列表被解释为int3_t（a、b、type 三个分量）结构体数组，逐键计算fbond = -2*k*dr/r，对两原子施加±delx*fbond等力，并调用Bond*重载的ev_tally_thr()累加能量与 virial。

归约逻辑（reduce_thr()的case THR_BOND分支）与模板化内核和 pair 情形几乎完全相同，这使得把任意 bonded 风格移植到 OpenMP 相对直接。angle/dihedral/improper 同理，参见 angle_cosine_omp.cpp、dihedral_charmm_omp.cpp 等；其中 CHARMM 角/二面体风格因含 1-3 成对项而使用组合标志THR_ANGLE|THR_CHARMM，reduce_thr()中为此设有专门分支，把eng_vdwl/virial_pair归约到全局 pair 对象上。

运行期配置：package omp 命令

编译只是前提。LAMMPS 输入脚本中还必须显式启用多线程，否则ThrOMP构造时会直接报错——fix_omp.cpp 中的检查：

fix = static_cast<FixOMP *>(lmp->modify->get_fix_by_id("package_omp")); if (!fix) lmp->error->all(FLERR, Error::NOLASTLINE, "The 'package omp' command is required for /omp styles");

对应的输入命令为：

package omp 4 neigh

第一个整数参数是每个 MPI 任务的线程数；传0表示"用omp_get_num_threads()探测到的全部线程"（源码中strcmp(arg[3],"0") == 0分支）。若与comm->nthreads不同，会调用omp_set_num_threads()重置并打印set N OpenMP thread(s) per MPI task。
可选关键字neigh yes/no控制邻居列表生成是否也多线程化（_neighbor标志）；未启用 OpenMP 编译时该命令只警告using 1 thread only并退回串行邻居列表。
package omp内部创建了ThrData数组、管理各 OMP 风格的注册顺序（last_omp_style/last_pair_hybrid字段用于决定何时执行最终的力归约，保证 hybrid 场景正确）。

使用方式即pair_style lj/cut/omp 2.5这类带/omp后缀的风格名。suffix_flag |= Suffix::OMP与注册名共同保证pair_style lj/cut（串行）与pair_style lj/cut/omp（线程）可以并存，用户按需在输入脚本中切换。

小结：编写新 OMP 风格的检查清单

综合文档与源码，为 LAMMPS OPENMP 包新增一个风格可按以下清单执行：

在src/OPENMP/新建xxx_omp.h/.cpp，类名如XxxOMP，双重继承Xxx（串行基类）与ThrOMP(lmp, THR_*)，构造中置suffix_flag |= Suffix::OMP，并在#ifdef PAIR_CLASS（或对应的BOND_CLASS等）分支中注册带后缀的风格名；
只 overridecompute()（bonded 风格还可顺带处理init()中需要线程化时机的部分）；respa_enable = 0显式关闭不支持的 RESPA；
compute()内按ev_init → omp parallel → loop_setup_thr → fix->get_thr(tid) → ev_setup_thr → eval<...>(ifrom, ito, thr) → reduce_thr的固定骨架编写，分区对象是邻居列表条目数或键列表长度；
把热循环抽成私有模板函数eval<EVFLAG, EFLAG, NEWTON_*>，用_noalias声明指针，力先写thr->get_f()，能量/virial 用对应类型重载的ev_tally_thr()；
若为多体势：改用 full 邻居列表，必要时用ev_tally3_thr()/ev_tally4_thr()等专用统计函数，在全局依赖点（归约、ghost 通信）用sync_threads()屏障并让 MPI 通信只落在#pragma omp master；
实现memory_usage()叠加memory_usage_thr()；同时确认 CMake（-DLAMMPS_OMP=yes或 preset）与传统Makefile.omp两种构建方式都能编译；
用package omp <nt> neigh输入脚本验证，并可结合compute pe/omp或逐原子量（compute atom/omp系列风格同样在本包内）核对线程结果与串行版一致。