ASE 径向分布与化学分波

本例使用集群现有 ASE 3.29.0、Python 3.11、NumPy 2.3.5 与 Matplotlib 3.11.2，处理已有真实 GPUMD/VASP 输出，不需要新提交 GPU 或 DFT 任务。下载完整结构、原始结果和可运行脚本，在新解压目录执行：

bash run.sh

analysis/ 保存 CSV、结果汇总和可继续使用的数据，figures/ 保存三张高清图。脚本只用一个 CPU 线程；重复运行请另解压一份，以免混入旧输出。

1. 先确认结构适合体相RDF#

这里读取GPUMD替位合金实际优化后的 512 原子 Si0.5Ge0.5 结构，Si、Ge 各 256 个，三个方向周期，体积 10892.421184 ų。它是一份经典势优化构型，不是有限温度平均，也不是 DFT 结构。

head -n 5 data/alloy.xyz

体相 RDF 以体积和数密度归一化。不要直接对带大量真空的二维晶胞、开放边界的异质界面用整盒体积，再把所得曲线当作体相 g(r)。本包专门使用充分大的周期体相合金，最大半径 8 Å，小于该胞的半盒高度。

2. 提取总RDF#

module load ase/3.29.0
python - <<'PY'
from ase.io import read
from ase.geometry.rdf import get_rdf
import numpy as np
a = read('data/alloy.xyz')
g, r = get_rdf(a, rmax=8.0, nbins=160)
np.savetxt('rdf-total.csv', np.c_[r, g],
           delimiter=',', header='r_A,g_r', comments='')
PY

rmax 是 Å，nbins=160 对应 0.05 Å 的箱宽，返回 r 为箱中心。当前版优先用 ase.geometry.rdf.get_rdf；旧 Analysis.get_rdf 已弃用，部分元素选择还存在已知问题，不建议继续照旧接口写新的化学分波程序。

3. 按元素对提取分波#

g_SiSi, r = get_rdf(a, 8.0, 160, elements=('Si', 'Si'))
g_SiGe, r = get_rdf(a, 8.0, 160, elements=('Si', 'Ge'))
g_GeSi, r = get_rdf(a, 8.0, 160, elements=('Ge', 'Si'))
g_GeGe, r = get_rdf(a, 8.0, 160, elements=('Ge', 'Ge'))

各分波按中心元素数和邻居元素数密度分别归一化。SiGe 与 GeSi 的归一化曲线应一致。对本例各占一半的成分：

g_total = 0.25 × (g_SiSi + g_SiGe + g_GeSi + g_GeGe)

不能直接把四个分波相加就称作总 RDF，也不能把分波的面积不经数密度换算直接当配位数。

4. 从积分得到配位数#

配位数为 ρ × ∫ g(r) 4πr²dr。脚本按每个球壳的精确体积计算积分，并用独立的全原子对最小镜像距离计数核对。积分到 3 Å，本例每个原子第一近邻均为 4，RDF 积分同样得到 4.000000。

head analysis/alloy-rdf.csv
cat analysis/RESULTS.json
ls figures

晶体单构型的峰很尖，改变箱宽会改变峰高，不能把峰高变化直接说成配位变了。应比较相同归一化下的峰位置与积分。对已有平衡轨迹，可用帧列表取平均：

g_average, r = get_rdf(frames, rmax=8.0, nbins=160)

用于真实热轨迹时，确认原子数、成分、体积与周期边界符合分析目的；变胞轨迹按每帧实际体积处理。进一步几何操作见ASE近邻与几何工具，原始合金模型见GPUMD合金优化。依据：ASE RDF实现与定义、旧Analysis接口说明。