OVITO晶体识别、化学键与结构因子
从理想晶体和苯分子出发，实际识别结构、计算键序和六元环，提取结构因子。
阅读本篇目录
下载输入、完整脚本和实际结果，在新目录解压后运行：
bash run.sh
分析使用OVITO3.16.1的Python接口，作图使用ASE3.29.0环境中的Matplotlib。共享环境提供预编译程序，用户无需自己编译。inputs/保存结构和轨迹，reference/保存参考数据，新结果写入work/。下面的单项命令在完整运行之后执行。
1. 先分清输入是什么#
本例的Al FCC、Fe BCC、Si金刚石都是明确晶格常数的理想几何；“冰”只构造立方冰的氧子晶格，不含氢，也没有进行成冰模拟。苯分子来自ASE几何数据库。它们适合判断分析工具的输入和输出，不能当成有限温度稳定性结论。
module purge
module load ovito/3.16.1
cd work
python examples.py classify
| 结构与原子数 | 方法 | 本例结果 |
|---|---|---|
| Fe BCC，128原子，a=2.87 Å | Ackland–Jones、VoroTop | 全部BCC |
| Al FCC，108原子，a=4.05 Å | CNA、PTM、VoroTop | 全部FCC |
| Si金刚石，216原子，a=5.43 Å | IdentifyDiamond | 全部立方金刚石 |
| 理想冰氧子晶格，216个O，a=6.35 Å | CHILL+，cutoff=3.5 Å | 全部立方冰型局域环境 |
| Al FCC | 中心对称参数，12近邻 | 最大约3.23×10⁻²⁹ Ų |
实际管线很短：
from ovito.io import import_file
from ovito.modifiers import PolyhedralTemplateMatchingModifier
p = import_file('inputs/fcc.xyz')
p.modifiers.append(PolyhedralTemplateMatchingModifier(output_orientation=True))
d = p.compute()
print(d.particles['Structure Type'])
print(dict(d.attributes))
更换算法时，输出数字编码也可能变化；按对应类的Type枚举解释，不能把所有方法的数字1一律理解为FCC。CHILL+分析水时应先只保留氧位置。VoroTop还需要filter，本包带上游BCC/FCC过滤器和许可；筛选器不是势函数。
2. 苯的成键、键序、键长与环#
python examples.py bonds
head results/benzene-bond-length.csv
head results/benzene-bond-angle.csv
from ovito.io import import_file
from ovito.modifiers import (CreateBondsModifier, BondOrderModifier,
FindRingsModifier, BondLengthDistributionModifier,
BondAngleDistributionModifier)
p = import_file('inputs/benzene.xyz')
p.modifiers.append(CreateBondsModifier(mode=CreateBondsModifier.Mode.CovalentRadius))
p.modifiers.append(BondOrderModifier(redistribute_bond_orders=True))
p.modifiers.append(FindRingsModifier(minimum_ring_size=3, maximum_ring_size=8))
p.modifiers.append(BondLengthDistributionModifier())
p.modifiers.append(BondAngleDistributionModifier())
d = p.compute()
print(d.particles.bonds['Bond Order'])
print(d.tables['ring-size-histogram'].xy())
得到12条键：6条C–H和6条C–C；芳香键重新分配后C–C键序为1.5，C–H为1。关闭redistribute_bond_orders时为整数单/双键表示。环统计只有1个六元环，键长直方图计数合计12，键角合计18。键序是基于几何与化学类型的推断，不能代替电子结构软件计算的键级。
元素名称必须是C、H等实际化学符号。若轨迹只有数字类型，先正确映射元素再使用共价半径建键。成键截断太大可能把非键合邻居也连起来。旧BondAnalysisModifier在本版本仍可计算，但新脚本分别使用键长与键角两个处理器更明确。
3. 从FCC位置得到结构因子#
python structure_factor.py
head results/fcc-direct-structure-factor.csv
from ovito.modifiers import StructureFactorModifier
p = import_file('inputs/fcc.xyz')
p.modifiers.append(StructureFactorModifier(
mode=StructureFactorModifier.Mode.Direct,
k_min=1, k_max=5, k_bins=200))
data = p.compute()
print(data.tables['structure-factor'].xy())
本例峰位为2.69、3.11、4.39 Å⁻¹，对应FCC的111、200、220反射。由a=4.05 Å计算的理论位置分别为2.68711、3.10281、4.38803 Å⁻¹；输出位置取直方图中心，存在0.02 Å⁻¹分箱精度。不要把本例归一化峰高直接解释为实验粉末衍射强度。
脚本还把前54个和后54个Al人为标记为A/B，只用于检查分项规则：总量等于AA+2AB+BB。最大数值差约4.44×10⁻¹⁶；这不是二元合金模型。对Na真实轨迹，python examples.py rdf还会生成Debye形式的结构因子；小体系、有限RDF截断和窗函数会影响曲线，不能仅凭峰的有无判断相变。
参考：OVITO结构因子、VoroTop过滤器。