OVITO属性计算与轨迹统计
计算、筛选和汇总粒子属性，冻结参考值，生成时间序列与轨迹线。
阅读本篇目录
下载输入、程序与实际结果，在新目录解压后先运行bash run.sh。新结果在work，参考数据在reference。分析环境为OVITO3.16.1，作图环境为ASE3.29.0。本文的单项命令在完整运行之后执行。
1. 先处理一帧的属性#
module purge
module load ovito/3.16.1
cd work
python examples.py properties
head results/R2-histogram.csv
示例Al FCC共108原子，a=4.05 Å、3×3×3常规胞。使用位置定义辅助属性R²=x²+y²+z²；它依赖原点，只演示如何生成和处理属性，不是均方位移。
from ovito.io import import_file
from ovito.modifiers import (ComputePropertyModifier,
HistogramModifier, ReducePropertyModifier, RemovePropertyModifier)
p = import_file('inputs/fcc.xyz')
p.modifiers.append(ComputePropertyModifier(
output_property='R2',
expressions=['Position.X^2+Position.Y^2+Position.Z^2']))
p.modifiers.append(HistogramModifier(property='R2', bin_count=20))
p.modifiers.append(ReducePropertyModifier(
property='R2', reduction_operation=ReducePropertyModifier.Operation.Mean))
d = p.compute()
print(d.tables['histogram[R2]'].xy())
print(dict(d.attributes))
p.modifiers.append(RemovePropertyModifier(properties=['R2']))
R²平均值112.7671875 Ų，直方图计数合计108。ReducePropertyModifier还可改为Sum、Min、Max、Median、Variance等操作；方差是样本方差，应用时注意归一化。RemovePropertyModifier删除整个属性，不能只删矢量的一个分量。
SpatialBinningModifier按空间区域汇总，本例沿x分6箱，逐箱x坐标求和再累加，与全体x坐标之和一致；不要把Sum直接叫作密度。Voronoi分析的108个原子体积相加为1793.613375 ų，与晶胞体积一致。
图左对一个FCC原子的12近邻建立配位多面体；图右是NumPy与Matplotlib从同一结构直接画出的属性散点图。本安装版的ScatterPlotModifier没有公开可配置x/y属性的Python字段，自动脚本使用这个等效绘图步骤，不依赖GUI状态文件。
2. 选择、扩展、删除与几何编辑#
python examples.py selection
python examples.py edit_geometry
python aliases.py
| 操作 | 实际用法 | 要检查的地方 |
|---|---|---|
| 按位置选原子 | ExpressionSelectionModifier(expression='Position.X < 6.075') |
本例选54/108原子 |
| 反选 | InvertSelectionModifier() |
两种选择数量之和108 |
| 删除所选 | DeleteSelectedModifier() |
删除的是Selection=1的原子 |
| 按类型选 | SelectTypeModifier(property='Particle Type', types={1}) |
类型1必须确实是目标元素 |
| 清空选择 | ClearSelectionModifier() |
本例直接移除Selection属性 |
| 扩展一个原子的近邻 | ExpandSelectionModifier(cutoff=3) |
FCC一层加中心共13原子 |
| 切片 | SliceModifier(normal=(1,0,0), distance=6.075, slab_width=3, select=True) |
默认选中薄片外侧的90原子；内侧18原子 |
| 复制 | ReplicateModifier(num_x=2, num_y=1, num_z=1) |
原子数、体积均翻倍 |
| 修改周期条件 | EditSimulationCellModifier(pbc=(False,False,False)) |
只改边界条件，坐标不变 |
| 合并与去重 | CombineDatasets→FindOverlappingParticles→DeleteSelected | 两份同一结构216原子，去重后108 |
注意区分粒子数组下标与Particle Identifier。ManualSelectionModifier插入时会捕获已有选择：本例先选ID=1，再改变上游表达式为前5个原子，捕获结果仍只有1个。自动化选指定ID通常直接用表达式；GUI交互切换不属于这里的独立Python流程。
3. 冻结值与两帧相减#
python examples.py time_analysis
head results/Na-MSD.csv
本例复用真实VASP Na轨迹：16原子、1000帧，每帧0.002 ps。先解包周期边界，再计算相对第0帧的位移，最后将位移平方在原子间取平均。
from ovito.modifiers import (UnwrapTrajectoriesModifier,
CalculateDisplacementsModifier, ComputePropertyModifier,
ReducePropertyModifier, TimeSeriesModifier)
p = import_file('inputs/Na-MD.xyz')
p.modifiers.append(UnwrapTrajectoriesModifier())
p.modifiers.append(CalculateDisplacementsModifier(minimum_image_convention=False))
p.modifiers.append(ComputePropertyModifier(
output_property='Displacement2',
expressions=['Displacement.X^2+Displacement.Y^2+Displacement.Z^2']))
p.modifiers.append(ReducePropertyModifier(property='Displacement2'))
key = next(k for k in p.compute().attributes if k.startswith('ReduceProperty'))
p.modifiers.append(TimeSeriesModifier(operate_on=key))
print(p.compute().tables['time-series'].xy())
第999帧的MSD为19.340255 Ų；脚本还用分数坐标逐帧解包，独立得到同一曲线。此处没有拟合扩散系数：16原子和约2 ps轨迹不足以支持普适的扩散结论。
FreezePropertyModifier(source_property='Position', destination_property='Initial Position', freeze_at=0)将初始值保留下来；DifferenceBetweenFramesModifier(operate_on='particles:particles/Position')只是坐标数值相减，不处理周期越界，需要物理位移时使用CalculateDisplacements。
4. 平滑、时间平均、轨迹线与合并轨迹#
同一time_analysis命令还执行以下步骤：
| 目的 | 设置 | 本例核对 |
|---|---|---|
| 平滑坐标 | 先Unwrap，再SmoothTrajectoryModifier(window_size=3, minimum_image_convention=False) |
第5帧等于第4/5/6帧坐标均值 |
| 平均RDF | RDF(cutoff=4, number_of_bins=80)→TimeAveraging(table:coordination-rdf) |
与逐帧RDF手工均值一致 |
| 画轨迹线 | 选择ID≤2→GenerateTrajectoryLines | 2个原子×1000帧，共2000个轨迹顶点 |
| 为拓扑加载轨迹 | LoadTrajectoryModifier的source.load | 第5帧坐标与原文件一致 |
平滑会改变短时振动信息，不能把平滑后的速度或位移直接代入原始动力学公式。时间平均前明确采样间隔、平衡段、原子ID和晶胞变化；本例保持一致原子ID和固定晶胞。
参考：OVITO Python API。