GPUMD轨迹的静态结构因子
用真实Si轨迹计算S(q)，提供ASE终端工具、数据与图片，并连接RDF、扩散、VACF、DOS和分组MSD。
阅读本篇目录
RDF看原子间距离，静态结构因子S(q)从倒空间观察原子排列。本篇用鸢算科技集群已有的7组真实GPUMD轨迹，配合ASE3.29.0和Python计算S(q)，不重复跑MD。原模拟为GPUMD5.9.1，Si_Fan_2019 Tersoff势，GCC13.3/CUDA12.9构建。
1. 下载、先看参考结果#
下载真实轨迹、工具、数据与图片，保存为gpumd-sq.zip：
unzip gpumd-sq.zip -d gpumd-sq
cd gpumd-sq
ls inputs reference-analysis reference-figures
head reference-analysis/source-summary.csv
包内完整XYZ保留原坐标和晶胞。reference-*是本集群实际后处理结果；重新计算写入独立results/analysis/figures。
2. 一条短命令处理一个轨迹#
在新的单例目录中：
module load ase/3.29.0
mkdir one
cd one
python ../sq.py ../inputs/07-structure/trajectory.xyz
head sq-result/sq.csv
换成自己的文件时，用python sq.py trajectory.xyz即可。输出第一列为含2π的物理波矢大小q，单位Å⁻¹；第二列为无量纲的原子数—原子数S(q)。配套脚本目前要求单元素、固定且与坐标轴对齐的立方周期晶胞，默认q上限3.8Å⁻¹。多组分、可变晶胞、非立方晶胞应使用相应算法，不能删掉检查后照搬。
完整重现7组对照时，在另一个新解压目录提交：
sinfo
sbatch job.slurm
模板按当前CPU分区56物理核整节点规则申请；在已分配的CPU作业内也可bash run.sh。此脚本只做CPU后处理，不调用GPUMD、不占GPU。单例与完整流程使用不同新目录，避免已有sq-result触发防覆盖检查。
3. 计算具体经过哪些步骤#
- ASE读取每帧的实际坐标、元素和晶胞。
- 按周期晶胞生成离散倒格矢，排除q=0作图。
- 每一帧先计算复数密度之和，再取模平方并除原子数。
- 对同一q长度的方向做壳平均，再对一条轨迹的帧平均。
- 同一目标温度两条独立初速度轨迹，先各自平均，再合并种子均值。
- 保存逐帧倒格矢数据、径向CSV、温度表与图片。
公式为S(q)=|sum_j exp(i q·r_j)|²/N。晶胞倒格矢包含2π，不同资料的q/k记号可能不同，比较横轴前先核对定义。不要先平均原子坐标再计算S(q)，那会丢掉轨迹中的结构信息。
4. 实际轨迹与原始输入#
| 原轨迹 | 原子数 | 保存帧数 | 生产段平均温度/K |
|---|---|---|---|
| 07-structure | 512 | 2 | 297.766 |
| liquid-si-3000-s51 | 1000 | 3 | 2963.259 |
| liquid-si-3000-s52 | 1000 | 3 | 3094.336 |
| liquid-si-3500-s51 | 1000 | 3 | 3534.077 |
| liquid-si-3500-s52 | 1000 | 3 | 3558.656 |
| liquid-si-4000-s51 | 1000 | 3 | 4059.353 |
| liquid-si-4000-s52 | 1000 | 3 | 4135.859 |
每个高目标温度有两个初速度种子，各只保留3帧；300K原轨迹保留2帧。目录里的liquid是原算例名，不能据此认定已经熔化；表中温度按整个生产段thermo.out计算，与目标温度及只统计后半段的均值不同。
产生4000K目标温度轨迹的原始run.in如下，仅供理解数据来源；本篇配套作业不重复执行它：
potential potential.txt
velocity 4000 seed 51
time_step 1
ensemble nvt_nhc 4000 4000 100
run 50000
ensemble nve
compute_msd 10 1000
compute_sdc 10 1000
compute_viscosity 10 1000
compute_rdf 8 160 500
dump_thermo 1000
dump_xyz 100000 trajectory.xyz
run 300000
5. 看图：晶体尖峰和高温轨迹#
左图比较512原子理想几何参考与实际约300K轨迹。理想几何只用于核对公式，不是另一次DFT或MD计算。对于该4×4×4金刚石常规晶胞超胞：
| 常规晶胞反射 | 超胞倒格矢整数 | 理想S(q) |
|---|---|---|
| q=0 | 0,0,0 | 512，作图排除 |
| 111 | 4,4,4 | 256 |
| 220 | 8,8,0 | 512 |
| 100 | 4,0,0 | 约0，消光 |
| 200 | 8,0,0 | 约0，消光 |
右图使用对数纵轴，3000/3500K目标温度轨迹仍出现明显晶体尖峰，4000K目标温度呈现更平缓的结构分布。这应与同批扩散/MSD结果和RDF合看，不能仅由一条曲线确定精确熔点。
图带是两个种子均值的范围，不是95%置信区间。帧数很少，512与1000原子的尺寸也不同；Bragg峰高度还受原子数和壳内方向平均影响，不能把峰高变化全部归因于温度。
6. 输出文件和进一步使用#
# 在完整工作流结束后的解压根目录
cat analysis/SUMMARY.json
head results/07-structure/radial-shells.csv
ls results/07-structure/per-frame.npz
CSV记录q壳、方向数和平均S(q)；NPZ保留倒格矢、逐帧结果和实际时间。程序还用独立的原子对余弦求和、整体平移和周期折回核对结果。
这里所有原子散射权重均为1，不是加入原子散射因子和仪器展宽后的实验XRD强度。q=0为N，不能直接拿这个值当热力学极限的压缩率。
7. 把结构与动力学方法衔接起来#
| 想看什么 | 原生入口或处理工具 | 实操资料 |
|---|---|---|
| RDF、键角分布 | compute_rdf、compute_adf | RDF与ADF |
| MSD、扩散与速度相关 | compute_msd、compute_sdc | 高温扩散与速度相关 |
| 振动态密度 | compute_dos | VDOS |
| 静态结构因子 | 本包ASE/Python后处理 | 本篇 |
| 多元素分组MSD | compute_msd的all_groups | LLZO分组实操 |
静态结构因子的定义可对照freud官方说明。本包实际使用ASE和配套Python脚本，没有调用freud库。正式研究需增加轨迹采样、独立种子与尺寸对照；本例用于展示从已有数据到结果的操作。