LAMMPS compute ke 命令完全指南:平动动能计算、单位换算与源码级原理剖析 科研科学计算高性能计算【免费下载链接】lammpsPublic development project of the LAMMPS MD software package项目地址https://gitcode.com/gh_mirrors/la/lammps点击查看免费下载本文系统讲解 LAMMPS 中compute ke命令的语法、物理含义与底层实现重点剖析它如何计算一组粒子的平动动能每个粒子贡献 ½mv²以及它与thermo_style中ke、etotal关键字由温度反推动能之间的微妙差异。同时结合仓库源码说明单位换算因子mvv2e的作用、MPI 并行归约流程、KOKKOS 加速变体以及与compute ke/atom、compute ke/rigid、compute erotate/sphere的搭配使用场景。读者读完可正确配置动能计算与热力学输出并理解其数值结果的物理依据。一、命令概览与语法compute ke是 LAMMPS 内置的计算命令之一用于计算指定组group内全部粒子的平动动能translational kinetic energy之和返回一个全局标量。命令定义在 compute_ke.rst 文档页支持 KOKKOS 加速变体ke/kk。compute ID group-ID ke其中ID计算对象的唯一标识符必须由字母、数字或下划线组成后续通过c_ID引用group-ID要计算动能的原子组只有组内原子参与求和ke本计算命令的风格名称。最简单的用法是compute 1 all ke表示对全部原子计算平动动能计算对象 ID 为1。更多关于compute命令通用语法ID、组、style 的约束可参考 compute 命令文档。参数校验从源码看该命令的参数个数是严格固定的。构造函数在 src/compute_ke.cpp 中ComputeKE::ComputeKE(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (narg ! 3) error-all(FLERR, Illegal compute ke command); scalar_flag 1; extscalar 1; }即参数必须恰好为 3 个ID、group-ID、ke多一个或少一个都会触发Illegal compute ke command错误。同时这里设置了两个关键标志scalar_flag 1本计算输出一个全局标量extscalar 1该标量是extensive广延量即随体系原子数增大而增大这与温度等 intensive强度量标量形成对比也决定了它在热力学输出中是否会被按原子数归一化详见下文输出信息一节。二、物理定义½mv²该计算定义的物理量是组内所有粒子的平动动能之和[ KE \sum_{i \in \text{group}} \frac{1}{2} m_i v_i^2 ]其中 (m_i) 与 (v_i) 分别是第 (i) 个粒子的质量与速度。注意这里的动能只包含平动部分不包含旋转动能对应compute erotate/sphere、振动或内部自由度等贡献。源码实现两条质量路径compute_scalar()在 src/compute_ke.cpp 中实现其核心逻辑如下double **v atom-v; double *rmass atom-rmass; double *mass atom-mass; int *mask atom-mask; int *type atom-type; int nlocal atom-nlocal; double ke 0.0; if (rmass) { for (int i 0; i nlocal; i) if (mask[i] groupbit) ke rmass[i] * (v[i][0]*v[i][0] v[i][1]*v[i][1] v[i][2]*v[i][2]); } else { for (int i 0; i nlocal; i) if (mask[i] groupbit) ke mass[type[i]] * (v[i][0]*v[i][0] v[i][1]*v[i][1] v[i][2]*v[i][2]); }可以看到两条并行路径per-atom 质量路径当体系使用mass命令为每个原子单独指定质量即atom-rmass非空时直接取rmass[i]按类型质量路径当使用默认的按原子类型质量mass数组时通过mass[type[i]]按原子类型查表获得质量。两者都累加 (m(v_x^2v_y^2v_z^2))注意此时尚未乘以 ½然后用mask[i] groupbit位掩码判断原子是否属于目标组只累加组内原子的贡献最后通过MPI_Allreduce将所有 MPI 进程的局部和归约得到全局总量MPI_Allreduce(ke, scalar, 1, MPI_DOUBLE, MPI_SUM, world); scalar * pfactor;其中的 ½ 因子和单位换算因子合并体现在pfactor中它在init()中初始化void ComputeKE::init() { pfactor 0.5 * force-mvv2e; }也就是说最终结果是0.5 * mvv2e * Σ m v²mvv2e负责把质量 × 速度²换算到当前单位制下的能量单位详见下一节。三、单位换算因子 mvv2e 的取值mvv2e是 LAMMPS 内部的质量-速度到能量换算因子其数值由 units 命令 选定的单位制决定在 src/update.cpp 的Update::set_units()中为各单位制赋值单位制mvv2e 值说明lj1.0所有量为约化单位无量纲real48.88821291² ≈ 2389.0能量单位 kcal/molmetal1.0364269e-4能量单位 eVsi1.0能量单位 Jcgs1.0能量单位 ergelectron1.06657236能量单位 Hartreemicro1.0能量单位 picoJpg·μm²/μs²nano1.0能量单位 10⁻²¹ J由此可以直观理解compute ke结果的单位在units metal下返回 eV在units real下返回 kcal/mol在units lj下返回约化能量单位。因此比较不同单位制下的动能输出时必须先统一单位。为什么需要 mvv2e在 LAMMPS 内部速度v以长度/时间存储质量m以质量单位存储两者的乘积m·v²并不天然等于能量单位。mvv2e正是把质量单位 ×长度/时间²转换为该单位制下能量单位如 eV、kcal/mol、J的换算系数。compute ke乘以0.5 * mvv2e后输出的标量才直接是标准能量单位下的平动动能。四、与 thermo_style 中 ke / etotal 关键字的本质区别这是理解compute ke最容易踩坑的地方。文档明确指出compute ke与热力学输出thermo_style中ke或etotal关键字给出的动能存在微妙差别。compute ke直接用公式 ½mv² 对组内原子逐粒子求和只含平动动能不含旋转等其它自由度热力学输出中的ke关键字由体系的温度反推动能即每个自由度贡献 ½k_B T 的能量因此它隐含了当前温度计算所用 compute 的自由度统计方式。这一逻辑在 src/thermo.cpp 的Thermo::compute_ke()中得到印证void Thermo::compute_ke() { dvalue temperature-scalar; dvalue * 0.5 * temperature-dof * force-boltz; if (normflag) dvalue / natoms; }即ke T × ½ × dof × k_B其中temperature-dof是温度计算中统计的自由度数。而compute temp的自由度定义为见 compute_temp 文档[ T \frac{2 E_{\mathrm{kin}}}{N_{\mathrm{DOF}} k_B}, \qquad N_{\mathrm{DOF}} n_{\mathrm{dim}} N_{\mathrm{atoms}} - n_{\mathrm{dim}} - N_{\mathrm{fix DOFs}} ]实现于 src/compute_temp.cppvoid ComputeTemp::dof_compute() { adjust_dof_fix(); dof domain-dimension * natoms_temp; dof - extra_dof fix_dof; if (dof 0.0) tfactor force-mvv2e / (dof * force-boltz); }其中extra_dof默认取维度数2 或 3对应周期性边界平移不变性去掉的自由度fix_dof是fix shake、fix rigid等约束 fix 冻结的自由度。什么情况下两者一致文档明确指出当温度由默认的compute temp即 LAMMPS 启动时自动创建的thermo_temp计算时compute ke的结果与热力学输出的ke相同。因为默认温度计算正好统计了体系全部平动自由度。但以下情况会导致二者分道扬镳换用其它温度计算方式如compute temp/partial扣除特定方向速度、compute temp/region仅统计某个区域、compute temp/com扣除质心平动、compute temp/deform、compute temp/profile等它们会扣除不同的非热速度成分、统计不同的自由度平动/转动从而与直接求和 ½mv² 的compute ke不一致存在约束 fix当体系包含fix shake、fix rigid等约束时compute temp会自动扣除被冻结的自由度而compute ke仍对组内所有原子的全部平动速度求和两者自然不同组范围不同热力学输出基于thermo_modify指定的温度 compute 的原子组而compute ke基于你自己指定的组。一句话总结想要纯粹的平动动能用compute ke想要与温度自洽的动能用热力学输出中的ke关键字。这也是文档将其列为Related commands并单独强调区别的原因。五、输出信息全局标量、extensive 属性与引用方式compute ke输出一个全局标量组内所有粒子的动能总和可以被任何需要全局标量输入的 LAMMPS 命令使用。常见的引用方式包括thermo_style custom ... c_ID把动能写入热力学输出注意此时不再归一化与否取决于thermo_modify norm的设置variable xxx equal c_ID赋值给变量用于条件判断或后续计算print、fix print、dump等输出命令。关于 LAMMPS 输出体系的整体介绍可参考 Howto output 文档。extensive 属性的含义文档明确该标量是extensive广延量单位是能量单位units 指定的单位制下的能量单位如 eV、kcal/mol、约化能量等。extensive 的影响主要体现在热力学输出中当thermo_style custom输出c_ID时若启用了thermo_modify norm yes默认在部分单位制下开启例如units lj下默认把所有广延量按原子数归一化动能会被除以原子数输出为每原子动能而温度这类 intensive 量不会被归一化。这一点在 units.rst 中也有说明thermo_style默认会把所有 extensive 量按原子数归一化行为可通过thermo_modify norm调整。因此当你把compute ke接入热力学输出并发现数值变小时应检查norm设置而不是怀疑计算错误。六、实际用法示例例 1把动能写入热力学输出units metal atom_style atomic lattice fcc 3.615 region box block 0 10 0 10 0 10 create_box 1 box create_atoms 1 box mass 1 63.546 velocity all create 300.0 4928459 compute 1 all ke thermo_style custom step temp c_1 etotal thermo 100 run 1000输出中c_1即为体系平动动能eVtemp为温度Ketotal为总能量势能 动能由温度反推的动能部分。例 2用变量追踪动能变化compute 1 all ke variable kineng equal c_1 fix 1 all print 100 ${kineng} file kinetic.out run 10000每 100 步将当前平动动能追加写入kinetic.out。例 3与转动动能组合有限尺寸粒子对atom_style sphere的有限尺寸粒子体系平动动能之外还需旋转动能此时配合使用compute 1 all ke compute 2 all erotate/sphere variable etot equal c_1c_2compute erotate/sphere计算球形粒子的旋转动能 ½Iω²I 为转动惯量ω 为角速度详见 compute_erotate_sphere.rst。注意该命令要求原子存储半径和角速度由atom_style sphere定义且组内粒子必须是球体或点粒子非球体粒子不允许。七、相关命令与扩展变体compute ke/atom逐原子动能如果需要每个原子的动能而不是总和使用compute ke/atom见 compute_ke_atom.rstcompute 1 all ke/atom它输出一个 per-atom 向量同样为能量单位组外原子的值为 0.0常配合dump custom输出空间动能分布或配合compute reduce等做二次汇总。compute ke/rigid刚体平动动能当体系用fix rigid/fix rigid/small把多个原子绑定为刚体时直接求和每个原子的 ½mv² 无法反映刚体的整体平动。RIGID 包提供compute ke/rigid见 src/RIGID/compute_ke_rigid.cpp它需要额外一个 fix-ID 参数fix 1 all rigid molecule compute 2 all ke/rigid 1其实现通过fix rigid的extract_ke()接口获得所有刚体质心平动动能之和见 src/RIGID/fix_rigid.cppdouble FixRigid::extract_ke() { double ke 0.0; for (int i 0; i nbody; i) ke masstotal[i] * (vcm[i][0]*vcm[i][0] vcm[i][1]*vcm[i][1] vcm[i][2]*vcm[i][2]); return 0.5*ke; }即 KE ½ Σ M_body V_cm²其中 M_body 是刚体总质量V_cm 是刚体质心速度。compute ke/rigid的init()中会校验给定的 fix 存在且风格以rigid开头否则报错。加速变体ke/kk 与 OpenMPcompute ke声明了 KOKKOS 加速变体ke/kk。在 src/KOKKOS/compute_ke_kokkos.cpp 中Kokkos 版本与串行版本逻辑完全一致同样的mask groupbit筛选、rmass/mass[type]双路径、MPI_Allreduce归约后乘pfactor只是将原子循环改写为Kokkos::parallel_reduce的 lambda 归约并借助atomKK-sync()确保速度、质量等数据在目标执行空间CPU/GPU上同步。因此在使用 KOKKOS 构建-k on时compute ke/kk可在 GPU 上并行归约但数值结果与串行版一致。此外compute ke/atom还提供/kk与/omp两个加速变体。八、使用注意事项与限制限制Restrictionscompute ke本身没有任何额外限制任何原子风格atom_style atomic、charge、sphere、molecular等都可直接使用速度必须已初始化动能计算依赖原子速度数组atom-v在使用前需通过velocity命令或fix langevin等初始化速度否则结果为零或未定义组范围决定求和范围compute ke只统计组内原子若后续用compute_modify改变组或使用dynamic选项配合原子增减需确认组定义符合预期与温度的自洽性问题如前所述只有默认compute temp且无约束 fix 扣自由度时compute ke才与热力学ke关键字数值一致否则二者不可混用、不可互相校验extensive 与归一化接入thermo_style custom时注意thermo_modify norm对广延量的归一化行为单位制输出单位为当前units单位制下的能量单位跨单位制对比结果时必须自行换算依据 src/update.cpp 中的mvv2e表。小结compute ke是 LAMMPS 中最基础也最常用的动能计算命令语法简单compute ID group-ID ke物理定义清晰平动动能 ½mv² 求和输出一个 extensive 的全局标量能量单位。它的实现src/compute_ke.cpp体现了 LAMMPS 计算的典型模式——组掩码筛选 MPI 归约 单位换算mvv2e并拥有 KOKKOS 加速变体。理解它与热力学输出ke关键字的区别、以及compute ke/atom、compute ke/rigid、compute erotate/sphere等扩展命令的适用场景是正确分析体系动能的起点。赞分享科研科学计算高性能计算【免费下载链接】lammpsPublic development project of the LAMMPS MD software package项目地址https://gitcode.com/gh_mirrors/la/lammps点击查看免费下载相关推荐LAMMPS 刚体旋转动能计算compute erotate/rigid 命令详解与源码剖析LAMMPS 刚体旋转动能计算compute erotate/rigid 命令详解与源码剖析 compute erotate/rigid 是 LAMMPS 中科研科学计算高性能计算LAMMPS compute erotate/sphere/atom 命令详解逐原子转动动能计算的用法与源码实现LAMMPS compute erotate/sphere/atom 命令详解逐原子转动动能计算的用法与源码实现 在 LAMMPS 中模拟含旋转自由度的球形颗科研科学计算高性能计算LAMMPS compute erotate/sphere球型粒子旋转动能的计算原理、实战用法与源码剖析LAMMPS compute erotate/sphere球型粒子旋转动能的计算原理、实战用法与源码剖析 本文围绕 LAMMPS 的 compute erot科研科学计算高性能计算上一篇Ryujinx终极Nintendo Switch模拟器完整指南与快速上手教程下一篇HACS实战指南高效解决Home Assistant社区商店的5大核心问题创作声明:本文部分内容由AI辅助生成(AIGC),仅供参考