LAMMPS 分组质心与转动统计

下载分组统计实例包，在已分配的兼容CPU资源中运行：

unzip lammps-chunk-observables-wave250.zip -d chunk-practice
cd chunk-practice
bash run.sh

使用ase/3.29.0提供Python绘图环境，计算模块为lammps/30Sep2026-8173m，CPU构建使用Intel oneAPI2024.2与Intel MPI2021.13。流程只用1个CPU线程，输出results-cpu/、analysis/与figures/。重新练习请用新解压目录，避免覆盖已有文件。

这是一个16个点粒子、4个分组的受控运动例子。采用LJ约化单位，粒子质量为1或2，分子ID为7、23、41、99；没有键约束、库仑相互作用，LJ势的epsilon特意设为0。电荷用于演示偶极统计。它帮助看懂命令、单位、数组与周期边界，不代表真实材料的相互作用模型，也不是扩散或温度平衡的生产计算。

1 chunk、group和molecule ID的关系#

group选择参与计算的原子，chunk/atom再把这些原子分成多个组。一个chunk可以是一条链、一个分子或一个空间网格。本例按输入中的molecule ID分组，并压缩非连续编号：

head -n 32 inputs/free_neutral/atoms.data
cat inputs/free_neutral/in.cpu

compute mols all chunk/atom molecule compress yes nchunk once ids once
compute prop all property/chunk mols id count

压缩后的行号是1–4，原分子ID仍是7、23、41、99，c_prop的两列分别保存原ID与原子数。每组有4个原子，总质量为6。

对追踪固定分子的MSD，分组身份须稳定。动态空间网格可以换原子，不能直接把某个网格行当作同一个扩散分子。

2 质心、速度与回转半径#

compute com all com/chunk mols
compute wrapped all com/chunk mols wrap yes
compute vcm all vcm/chunk mols
compute gyr all gyration/chunk mols
compute gt all gyration/chunk mols tensor
compute shape all gyration/shape/chunk gt

质心与回转半径都按质量加权。本例原子质量不同，不能直接用四个坐标的算术平均代替质心。

head -n 14 results-cpu/free_neutral/prop.dat
head -n 14 results-cpu/free_neutral/com.dat
head -n 14 results-cpu/free_neutral/gyr.dat

初始每组Rg=0.67288764。第7号分子在约化时间1时变为0.70331556：粒子自由运动且没有刚性约束，形状会改变；molecule ID本身不会自动让分子保持刚体。

3 跨周期边界时读取unwrapped坐标#

23号分组初始质心x=49.8，质心速度vx=2，盒长为50；时间1时连续坐标为51.8，而包回盒内的坐标为1.8。两者表示同一个位置。

head -n 25 results-cpu/free_neutral/trajectory.dump
head -n 14 results-cpu/free_neutral/wrapped.dat

xu/yu/zu结合image flags保持跨边界分子的连续几何。本例输入中已经有部分原子跨过x=50，因此不能为图省事把所有image flags清零，否则Rg、惯性和偶极都可能被盒长污染。

图中连线仅帮助辨认4个点，不表示输入存在化学键。

4 惯性、角动量、等效角速度与力矩#

compute inertia all inertia/chunk mols
compute ang all angmom/chunk mols
compute omega all omega/chunk mols
compute torque all torque/chunk mols

angmom/chunk计算围绕各组质心的角动量，omega/chunk通过Iω=L求等效角速度。当前例子用点粒子；有限尺寸粒子的自转贡献应另按所用版本和atom style读取。

注意两种张量的列序不同：gyration/chunk tensor为xx、yy、zz、xy、xz、yz；inertia/chunk为xx、yy、zz、xy、yz、xz。不要复制一套列名用于所有张量。

head -n 9 results-cpu/free_neutral/inertia.dat
head -n 9 results-cpu/free_neutral/ang.dat
head -n 9 results-cpu/forced_neutral/torque.dat

初始第7号分组惯性对角元为1.72833333、1.70166667、2.00333333。恒力例子对每个原子施加(0.02,−0.03,0.01)，各原子质量不同，不能把它误读为相同加速度；由相对质心位置叉乘原子力得到组力矩。

5 把全局数组写成可读文件#

fix out_com all ave/time 10 1 10 c_com[*] file com.dat mode vector format " %.16g"
fix out_gyr all ave/time 10 1 10 c_gyr file gyr.dat mode vector format " %.16g"

每10步输出一次，每次只采1点，所以这里没有时间平滑。com/chunk是多列数组，用c_com[*]；标量回转半径按chunk组成全局向量，用c_gyr。自定义数字格式前保留空格，才能分开相邻数值列。

形状教程的三个特征值在本例输出为从大到小排列；后面三列为非球形度、非圆柱度和相对形状各向异性。别把特征值当作坐标长度，它们量纲是长度平方。

更多数据流操作见分组温度、MSD与数据传递，定义参考LAMMPS官方chunk说明。