OVITO-Python 处理 LAMMPS 轨迹：统计原子 X 方向密度分布及扩展

注意：ProbabilityDensity 是'概率密度'，不是物理意义上的'数密度（#/体积）'。要得到真正数密度，仍建议用 $V_{\text{bin}}$ 做体积归一化。

三、环境准备：安装 OVITO Python 模块

两种常见方式：

方式 A：pip（普通 Python 环境）

PyPI 页面列出了支持的平台与 Python 版本范围（例如 Windows/Linux x86_64、Python 3.10–3.13 等）。

方式 B：conda（推荐给 conda 用户）

OVITO 官方给了 conda 安装命令（含 ovito 模块；并提示 conda 包也可能带桌面端组件）。

四、数据准备：你能统计什么，取决于 dump 里有什么

4.1 统计密度（位置分布）最低要求

dump/xyz 至少要有：

• Particle Type（或 type）
• Position.X/Y/Z

如果你用的是 legacy XYZ（列含义不固定），需要用 import_file(..., columns=[...]) 显式告诉 OVITO 每一列是什么含义，例如：

pipeline = import_file("dump.xyz", columns=["Particle Type", "Position.X", "Position.Y", "Position.Z"])

OVITO 的 import_file() 会自动识别格式并创建 pipeline。

4.2 想统计速度/势能/应力：必须提前写进 dump

• 速度：dump 里要有 vx vy vz
• 势能：需要 LAMMPS compute pe/atom 并 dump custom 输出该列。
• 应力：需要 compute stress/atom（或相关命令）并输出到 dump；LAMMPS 文档也给过'用 per-atom stress 做 1D 压力 profile'的示例思路。

关键提醒：后处理不能凭空'算出'你没输出的 per-atom 势能/应力。OVITO 能做的是：对'文件里已有的每原子属性'进行统计、分箱、平均。

五、最小可用：Type=1 原子沿 x 的直方图（增强版）

原始代码非常接近正确，但建议补上两点：

• 显式开启 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]")

• SelectTypeModifier 的用法：按某个 property（例如 Particle Type）选择指定类型。
• HistogramModifier 的关键参数：property、bin_count、only_selected、fix_xrange/xrange_start/xrange_end、normalization。

'key=...到底是什么？'

很多 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",
)

八、常见坑（新手必看）

1. XYZ 列映射没写：legacy XYZ 必须用 columns=[...] 指定列含义，否则 OVITO 不知道哪列是 type、哪列是 x/y/z。
2. 设置了 xrange_start/end 但没生效：要 fix_xrange=True。
3. 把'计数'当成'密度'：真正数密度要除以 bin 体积（或至少除以 bin 宽度做 1D 线密度）。
4. 想统计势能/应力但 dump 没有：后处理救不了；只能回到 LAMMPS 增加输出（如 compute pe/atom、compute stress/atom）再跑或重 dump。
5. Key 写错：先 print(data.tables.keys()) 看看实际表名，再填 export_file(..., key=...)。OVITO 文档给过这种'按 key 导出表'的范式示例。