Python算法
OVITO-Python 处理 LAMMPS 轨迹：统计原子 X 方向密度分布及扩展
介绍如何使用 OVITO 的 Python 模块对 LAMMPS 模拟后的轨迹文件进行后处理。通过导入数据、选择特定类型原子、使用直方图修饰器（HistogramModifier）沿指定方向分箱统计，可计算数密度、速度或势能分布。文章涵盖了环境安装、数据准备、单帧统计、多帧平均以及常见错误排查，并提供了完整的代码示例。
介绍如何使用 OVITO 的 Python 模块对 LAMMPS 模拟后的轨迹文件进行后处理。通过导入数据、选择特定类型原子、使用直方图修饰器（HistogramModifier）沿指定方向分箱统计，可计算数密度、速度或势能分布。文章涵盖了环境安装、数据准备、单帧统计、多帧平均以及常见错误排查，并提供了完整的代码示例。
在分子动力学（MD）模拟的后处理中，沿某个方向的空间分布（profile）是非常常见的分析需求，例如数密度/质量密度、速度分布、应力（压力）分布、电荷分布、势能分布等。理想情况下，我们可以在 LAMMPS 的 in 文件里用 compute chunk/atom + fix ave/chunk 在运行过程中直接对空间分箱（binning）并做时间平均；而且 fix ave/chunk 还能对多种 per-atom 量（如速度、力、电荷、势能、应力等）做按 chunk 的汇总/平均输出。
但很多时候模拟已经跑完，脚本里没写这些 profile 计算，手头只剩下轨迹文件（dump/xyz/lammpstrj）。这时再回去重跑既费时又不必要——更合适的做法是：使用后处理工具从轨迹中重新统计分布。
OVITO 的 Python 模块提供了一套非常清晰的'数据流水线（pipeline）'机制：导入数据（import_file）→ 依次应用 modifiers → compute() 得到结果 → export_file 导出，非常适合把后处理流程脚本化、自动化。
在 LAMMPS 里，密度/电荷/速度/应力等空间分布，本来可以在 in 文件里通过 chunk + 平均直接算出来（典型做法是 compute chunk/atom + fix ave/chunk），并输出为 profile 文件。LAMMPS 官方文档也明确：chunks 可以按空间 bin 切分，并对每个 chunk 统计包括密度、速度、应力、势能等多种量。
但现实经常是：
计算已经跑完了；当时没有在
in文件里写 profile 的计算；现在只剩 dump/xyz 轨迹文件。
这时就需要后处理。OVITO 的 Python 模块非常适合做这一类'按空间分箱统计'的任务：读取轨迹 → 选择原子子集 → 做直方图/分箱 → 导出表格。
我们把盒子在 x 方向切成 bin_count 个小区间（bin），对每个 bin：
其中 $L_y, L_z$ 是盒子在 y,z 方向长度（假设截面积不随时间变化或你逐帧计算）。
OVITO 的 HistogramModifier 支持多种归一化方式：
注意：
ProbabilityDensity是'概率密度'，不是物理意义上的'数密度（#/体积）'。要得到真正数密度，仍建议用 $V_{\text{bin}}$ 做体积归一化。
两种常见方式：
PyPI 页面列出了支持的平台与 Python 版本范围（例如 Windows/Linux x86_64、Python 3.10–3.13 等）。
OVITO 官方给了 conda 安装命令（含 ovito 模块；并提示 conda 包也可能带桌面端组件）。
dump/xyz 至少要有：
如果你用的是 legacy XYZ（列含义不固定），需要用 import_file(..., columns=[...]) 显式告诉 OVITO 每一列是什么含义，例如：
pipeline = import_file("dump.xyz", columns=["Particle Type", "Position.X", "Position.Y", "Position.Z"])
OVITO 的 import_file() 会自动识别格式并创建 pipeline。
vx vy vzcompute pe/atom 并 dump custom 输出该列。compute stress/atom（或相关命令）并输出到 dump；LAMMPS 文档也给过'用 per-atom stress 做 1D 压力 profile'的示例思路。关键提醒：后处理不能凭空'算出'你没输出的 per-atom 势能/应力。OVITO 能做的是：对'文件里已有的每原子属性'进行统计、分箱、平均。
原始代码非常接近正确，但建议补上两点：
fix_xrange=True，否则 xrange_start/end 可能不生效。normalization（默认是绝对计数）。from ovito.modifiers import HistogramModifier, SelectTypeModifier
from ovito.io import import_file, export_file
pipeline = import_file("dump.xyz") # 或 dump.lammpstrj 等
# 选出 Type=1
pipeline.modifiers.append(
SelectTypeModifier(operate_on="particles", property="Particle Type", types={1})
)
# 沿 x 分箱做直方图
hist = HistogramModifier(
operate_on="particles",
property="Position.X",
bin_count=15,
only_selected=True,
)
hist.fix_xrange = True
hist.xrange_start = 0
hist.xrange_end = 50
# 绝对计数（默认就是这个）
hist.normalization = HistogramModifier.NormalizationMode.AbsoluteCount
pipeline.modifiers.append(hist)
# HistogramModifier 会生成一个 DataTable，名字通常类似 histogram[Position.X]
export_file(pipeline, "density.txt", "txt/table", key="histogram[Position.X]")
很多 modifier 会把结果以 DataTable 的形式放进 pipeline 输出里，然后 export_file(..., key=...) 指定导出哪张表。OVITO 文档示例明确展示了这种模式（比如导出 key='binning' 的表）。
实战建议：你可以先打印一下有哪些表：
data = pipeline.compute()
print(list(data.tables.keys()))
然后把 key 改成你实际看到的名字。
上面的 density.txt 本质是'每个 bin 的计数'。要变成数密度（#/体积），你需要盒子尺寸 $L_y, L_z$（以及 bin 宽度 $\Delta x$）。
下面给一个新手友好的脚本：导出两列：bin_center 与 number_density
import numpy as np
from ovito.modifiers import HistogramModifier, SelectTypeModifier
from ovito.io import import_file
pipeline = import_file("dump.xyz")
pipeline.modifiers.append(SelectTypeModifier(property="Particle Type", types={1}))
hist = HistogramModifier(property="Position.X", bin_count=15, only_selected=True)
hist.fix_xrange = True
hist.xrange_start = 0
hist.xrange_end = 50
hist.normalization = HistogramModifier.NormalizationMode.AbsoluteCount
pipeline.modifiers.append(hist)
data = pipeline.compute()
# 取出 histogram 表
table = data.tables["histogram[Position.X]"]
# 不同版本表列名可能略不同，最稳妥是打印 table.columns
# print(table.columns)
# 这里用'通用思路'：第一列通常是 bin center，第二列是 count
bin_center = np.array(table[:, 0])
count = np.array(table[:, 1])
# 盒子尺寸：data.cell 是模拟盒信息（假设正交盒）
# 若是三斜晶胞，需要更谨慎地取横截面积（这里先给初学者版）
Ly = data.cell[:, 1].length
Lz = data.cell[:, 2].length
dx = (hist.xrange_end - hist.xrange_start) / hist.bin_count
Vbin = dx * Ly * Lz
number_density = count / Vbin
np.savetxt(
"number_density_x.txt",
np.column_stack([bin_center, number_density]),
header="bin_center number_density(#/volume)",
)
如果你是 NPT、盒子会变化：建议逐帧计算并平均（下一节给模板）。
OVITO 的轨迹是'按需加载'的：不会一次读入所有帧，而是你
compute(frame)时读那一帧。 因此你可以这样做'逐帧累加 + 平均'：
import numpy as np
from ovito.io import import_file
from ovito.modifiers import SelectTypeModifier, HistogramModifier
pipeline = import_file("dump.xyz")
pipeline.modifiers.append(SelectTypeModifier(property="Particle Type", types={1}))
hist = HistogramModifier(property="Position.X", bin_count=100, only_selected=True)
hist.fix_xrange = True
hist.xrange_start = 0
hist.xrange_end = 50
pipeline.modifiers.append(hist)
sum_count = None
bin_center = None
nframes = pipeline.num_frames
for frame in range(nframes):
data = pipeline.compute(frame)
table = data.tables["histogram[Position.X]"]
bc = np.array(table[:, 0])
ct = np.array(table[:, 1])
if sum_count is None:
sum_count = ct.copy()
bin_center = bc
else:
sum_count += ct
avg_count = sum_count / nframes
np.savetxt(
"avg_count_x.txt",
np.column_stack([bin_center, avg_count]),
header="bin_center avg_count",
)
columns=[...] 指定列含义，否则 OVITO 不知道哪列是 type、哪列是 x/y/z。fix_xrange=True。compute pe/atom、compute stress/atom）再跑或重 dump。print(data.tables.keys()) 看看实际表名，再填 export_file(..., key=...)。OVITO 文档给过这种'按 key 导出表'的范式示例。
微信公众号「极客日志」，在微信中扫描左侧二维码关注。展示文案：极客日志 zeeklog
使用加密算法（如AES、TripleDES、Rabbit或RC4）加密和解密文本明文。 在线工具，加密/解密文本在线工具，online
解析常见 curl 参数并生成 fetch、axios、PHP curl 或 Python requests 示例代码。 在线工具，curl 转代码在线工具，online
将字符串编码和解码为其 Base64 格式表示形式即可。 在线工具，Base64 字符串编码/解码在线工具，online
将字符串、文件或图像转换为其 Base64 表示形式。 在线工具，Base64 文件转换器在线工具，online
将 Markdown（GFM）转为 HTML 片段，浏览器内 marked 解析；与 HTML转Markdown 互为补充。 在线工具，Markdown转HTML在线工具，online
将 HTML 片段转为 GitHub Flavored Markdown，支持标题、列表、链接、代码块与表格等；浏览器内处理，可链接预填。 在线工具，HTML转Markdown在线工具，online