ABACUS 轨迹的径向分布与结构因子

下载CANDELA源码、完整INPUT、实际ABACUS轨迹、原生结果和绘图脚本，解压进入附件根目录。先提取已有结果：

bash analyze.sh

图在 reference/figures，对照数据在 reference/analysis。附件里的两段Al32轨迹来自鸢算科技集群真实ABACUS3.10.1计算，各81帧、步长0.5 fs，共40 fs；一段NVE、一段Nose–Hoover链控温。它们用于教学，保持了晶体结构，不能当作已经平衡的液态或长时间扩散样本。

1 使用哪个程序，怎样重新运行#

CANDELA是轨迹分析程序，读取ABACUS的 MD_dump，不重新计算DFT。这里使用官方源码固定提交 9d4befd83c66b23463b8bac8939b98bb1516f72e，GNU C++/Fortran13.3，通用x86-64、串行、不开CUDA/OpenMP。附件带源码和许可，在自己的算例目录编译，不需要改动公共软件。

新解压目录在已分配的计算资源中运行：

bash run.sh

也可以提交附件模板：

sbatch submit.slurm

模板选P100主机分配资源，但本例只执行单线程CPU后处理，不是GPU加速测试。它依次编译CANDELA、建立独立 work、运行6个分析条件，再提取和画图。已存在 work 或 candela-build 时请新解压，避免覆盖结果。

希望逐步执行，可以在已分配资源中运行下面这些命令；它与 run.sh 二选一：

bash build.sh
python3 -B setup.py
python3 -B run_cases.py

CANDELA的输入也叫 INPUT，不要在ABACUS原始计算目录直接运行它。 本附件把轨迹、分析输入、参考结果、新结果分开保存。

2 计算径向分布函数#

查看实际PDF输入：

cat inputs/al-nve-pdf/INPUT

calculation pdf
system Al
geo_in_type ABACUS
geo_directory ../../trajectories/al-nve.MD_dump
geo_1 0
geo_2 80
geo_interval 1
geo_ignore 0
ntype 1
natom 32
natom1 32
id1 Al
id2 Al
msd_dt 0.0005
geo_out pdf.txt
rcut 3.9
dr 0.05
struf_dgx 0.05
struf_ng 100

geo_in_type ABACUS读取原生格式，geo_directory指向复制的轨迹；0到80包含81帧。rcut=3.9 Å小于本例8.1 Å立方盒的一半，dr=0.05 Å决定直方图宽度。msd_dt在这里一并写明为0.0005 ps，对MSD尤其重要。换轨迹时先检查元素数量、盒长、输出帧间隔和索引。

建立work并完成编译后，单独执行一个分析目录：

cd work/cases/al-nve-pdf
../../../candela-build/bin/candela

上述手动方式直接写当前目录输出；若使用 run_cases.py 已完成同一目录，不要重复运行。完整流程建议用 run.sh，它会同时记录执行信息并运行两个轨迹的三种分析。

head -n 8 reference/cases/al-nve-pdf/pdf.txt

pdf.txt第1列是壳层中心r/Å，第2列是g(r)。本版对单元素使用N(N−1)有限原子数归一化，32个原子时要保留这项修正。第3列是程序的径向积分累计量，没有乘数密度，不能直接命名为配位数；若要配位数，须明确密度、有限N约定与积分范围后再计算。

图中实线为原生CANDELA，点线由实际坐标独立做周期最近像、成对计数和球壳归一化。NVE最大g(r)差约4.43×10⁻⁵，NHC约1.48×10⁻⁵，主要受原生文本输出精度影响。短轨迹展宽后的峰仍体现晶体局部结构，不是液体径向分布的收敛判断。

3 直接计算静态结构因子#

cat inputs/al-nve-ssf/INPUT

calculation ssf
system Al
geo_in_type ABACUS
geo_directory ../../trajectories/al-nve.MD_dump
geo_1 0
geo_2 80
geo_interval 1
geo_ignore 0
ntype 1
natom 32
natom1 32
id1 Al
id2 Al
msd_dt 0.0005
struf_dgx 0.775702198941453
struf_dgy 0.775702198941453
struf_dgz 0.775702198941453
struf_ng 3
ssf_out ssf.txt

本例用盒子倒格矢间距2π/L，在三个方向各取−3到3的整数网格；原生程序按|q|相同的壳层平均。定义是 S(q)=〈|Σ exp(iq·rj)|²/N〉，q为Å⁻¹。0到80每帧都读取，先保持 geo_interval 1，避免跳帧设置与原生顺序读取逻辑不匹配。

head -n 8 reference/cases/al-nve-ssf/ssf.txt

两个轨迹的原生S(q)与独立复数密度振幅求和最大差约2.9×10⁻⁵。原生内部浮点累加及打印精度会留下微小误差。这里比较的是同一实际倒格矢集合和壳层，不能拿另一种盒子、不同q采样直接对齐。

PDF目录也会生成 ssf.txt，那是由有限半径g(r)变换得到的另一条曲线，受截断和积分影响；本节图和表读取的是 al-nve-ssf/ssf.txt，即坐标直接计算的结果，不把两个同名输出混在一起。

4 重新提取与画图#

module load ase/3.29.0
python -B analyze.py reference
python -B plot_results.py reference

处理你自己的新结果时把 reference 换成 work。更换原子数、元素或变胞轨迹不能只改路径：本附件独立核对脚本针对Al32固定正交晶胞，自己的体系应修改相应检查和归一化。更多方法见CANDELA官方文档与固定源码仓库。